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Bottomonium spectrum at order v^ from domain-wall lattice QCD: 
precise results for hyperfine splittings 
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The bottomonium spectrum is computed in dynamical 2+1 flavor lattice QCD, using NRQCD 
for the b quarks. The main calculations in this work are based on gauge field ensembles generated 
by the RBC and UKQCD collaborations with the Iwasaki action for the gluons and a domain-wall 
action for the sea quarks. Lattice spacing values of approximately 0.08 fm and 0.11 fm are used, 
and simultaneous chiral extrapolations to the physical pion mass are performed. As a test for gluon 
discretization errors, the calculations are repeated on two ensembles generated by the MILC col- 
laboration with the Liischer-Weisz gauge action. Gluon discretization errors are also studied in a 
lattice potential model using perturbation theory for four different gauge actions. The nonperturba- 
tive lattice QCD results for the radial and orbital bottomonium energy splittings obtained from the 
RBC/UKQCD ensembles are found to be in excellent agreement with experiment. To get accurate 
results for spin splittings, the spin-dependent order-u® terms are included in the NRQCD action, 
and suitable ratios are calculated such that most of the unknown radiative corrections cancel. The 
cancellation of radiative corrections is verified explicitly by repeating the calculations with different 
values of the couplings in the NRQCD action. Using the lattice ratios of the 5-wave hyperfine and 
the IP tensor splitting, and the experimental result for the IP tensor splitting, the 15* hyperfine 
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PACS numbers: 12.38. Gc, 14.40.Pq 



I. INTRODUCTION 



The low-lying radial and orbital energy splittings in bottomonium are well known from experiment and well under- 
stood theoretically. Spin-dependent energy splittings pose a significantly greater challenge. On the experimental side, 
the observation of the S = states is very difficult. So far, only the r]b{lS) has been found. The weighted average of 
Cn the results from [iHS] gives a value of 69.3 ± 2.9 MeV for the T(IS') — r]i,{lS) hyperfine splitting. 

^ Calculations of this hyperfine splitting using perturbative QCD gave significantly lower values around 40 MeV [IHS] . 

Mixing of the 775 with a light CP-odd Higgs boson has been suggested as a possible explanation of this discrepancy [7] . 
1^^ To definitely answer the question whether QCD alone is able to correctly predict the T(IS') — rjii{lS) mass difference, 
^y-^ precise non-perturbative calculations from lattice QCD are required. 

• Presently, lattice calculations with dynamical light quarks are performed at lattice spacings that are too coarse to 

U-K resolve the Compton wave length of the b quark, and therefore special heavy-quark techniques are required. Two 

f-~^ such techniques are the Fermilab method [S] and lattice nonrelativistic QCD (NRQCD) [51 [10]. The bottomonium 
y—{ spectrum has also been calculated using relativistic actions, on anisotropic lattices |11) and on very fine, small lattices 
|12j . but so far without dynamical light quarks. 

With the Fermilab method, the heavy quark is implemented by an improved Wilson-like action, where the param- 
eters are tuned such that heavy-quark discretization errors are reduced. The Fermilab method has the advantage 
^ over NRQCD that continuum extrapolations can be performed safely. In the simplest case, only the mass parameter 
in the action is adjusted such that the kinetic mass of a heavy meson agrees with experiment. This method, in 
combination with MILC gauge field configurations generated with the Liischer-Weisz gluon action and 2-1-1 flavors of 
rooted staggered sea quarks |13j . has been employed in |14j to calculate the bottomonium and charmonium spectra. 
In the continuum limit, the T(15') — r]i,{lS) splitting was found to be 54 ± 12 MeV. 

A version of the Fermilab action with three tuned parameters was used in Ref. |15j to calculate bottomonium 
masses at one lattice spacing using gauge field ensembles generated by the RBC and UKQCD collaborations using 
the Iwasaki gluon action and 2-1-1 flavors of domain-wall sea quarks jl6) . 

In contrast to the Fermilab method, lattice NRQCD is based on the direct discretization of an effective field theory 
for heavy quarks, in which an expansion in powers of the heavy-quark velocity v is performed [SlIlOj. For bottomonium, 
one has f^ ss 0.1 [T7]. With NRQCD, it is required that the UV cut-off provided by the lattice is lower than the 
heavy-quark mass: one must have ami, ^ 1, where a is the lattice spacing. Discretization errors can be removed 
through Symanzik improvement. The bottomonium spectrum was calculated using improved lattice NRQCD of order 
-y'* on the MILC gauge field ensembles in [TB]. Good agreement with experiment was seen for the radial and orbital 
energy splittings. However, the results for the spin-dependent energy splittings suffered from large uncertainties due to 



missing radiative and higher-order relativistic corrections in the NRQCD action used there. Spin-dependent sphttings 
are an effect of order u^, and therefore v^ corrections are significant. These v^ corrections were included in earher 
calculations [TM2^ and found to reduce the bottomonium spin splittings by 10-30%. Radiative corrections to the 
spin-dependent couplings in the NRQCD action are expected to be of order as (~20-30%), and are still unknown. 
However, as will be demonstrated in this work, these radiative corrections largely cancel in suitable ratios of spin 
splittings. 

In the following, a new calculation of the bottomonium spectrum in lattice QCD with 2+1 flavors of dynamical sea 
quarks is presented. For the b quarks, an improved NRQCD action including the spin-dependent order-w^ terms is used 
(for comparison, results obtained without these terms are also shown) . By varying the couplings of the leading-order 
spin-dependent terms in the action, the cancellation of radiative corrections is demonstrated nonperturbatively for 
the ratio of the 5'-wave hyperfine and IP tensor splittings, as well as the ratio of the 25* and IS* hyperfine splittings. 
For these quantities, results with unprecedented precision are obtained here. 

The main calculations in this paper are done on RBC/UKQCD gauge field ensembles, which were generated with 
the Iwasaki gluon action and a domain-wall sea quark action. The calculations are performed at lattice spacing values 
of approximately 0.08 fm (with spatial lattice size L = 32) and 0.11 fm (for both L = 24 and L = 16). This work is 
an extension of the first calculation by the author that was using only the L = 24 ensembles and only the w^ action 
[23j . By including the finer L ^ 32 ensembles, discretization errors can now be studied directly, and by including 
the L = 16 ensembles with their smaller box size, the size of finite-volume effects can be estimated. The data at 
L = 24 are also reanalyzed with improved methods leading to smaller statistical errors. With the better accuracy, 
chiral extrapolations to the physical pion mass are now possible. 

Systematic errors caused by the lattice NRQCD action are estimated using power counting. In order to study 
discretization errors caused by the lattice gluon action, two approaches are used here. First, the radial and orbital 
energy splittings are calculated in a lattice potential model, using the static quark-antiquark potential derived from 
the gluon action in lattice perturbation theory. A comprehensive study of the scaling behavior is presented for four 
different gluon actions (Plaquette, Liischer-Weisz, Iwasaki, and DBW2). This model is however limited to the tree 
level, and to radial and orbital energy splittings only. In order to go beyond tree-level and include spin splittings, 
the nonperturbative lattice QCD calculations of the bottomonium spectrum are repeated on gauge field ensembles 
generated with the Liischer-Weisz gluon action by the MILC collaboration, and a detailed comparison to the results 
from the RBC/UKQCD ensembles (with the Iwasaki action) is made. 

This paper is organized as follows: the lattice methods and parameters are described in Sec. |TTj The "speed of 
light" is studied in Sec. |III A| In Sec. |IIIB[ the results for the radial and orbital energy splittings are presented. 



followed by the spin-dependent energy splittings in Sec. Ill C The conclusions are given in Sec. |IV| A simple analysis 



of autocorrelations is described in Appendix |Aj the tuning of the 6-quark mass is discussed in Appendix [Bj and 
tables with results in lattice units can be found in Appendix O The lattice potential model calculations of gluon 
discretization errors and the comparison of nonperturbative results from the MILC and RBC/UKQCD ensembles are 
presented in Appendix [P] 

II. METHODS 
A. Lattice actions and parameters 

The calculations in this work are based on gauge field ensembles that include the effects of dynamical up- down- 
and strange sea quarks (with ?7i„ = rrid, in the following denoted as mi). The ensembles used for the main part of 
the calculations were generated by the RBC and UKQCD collaborations [TBI [24]. The sea quarks are implemented 
with a domain wall action |25H27j , which yields an exact chiral symmetry when the extent L5 of the auxiliary fifth 
dimension is taken to infinity. The gluons are implemented with the Iwasaki action 28-30 , which suppresses the 
residual chiral symmetry breaking at finite L5 |31j . The form of the Iwasaki action can be found in Appendix pi 
where the discretization errors associated with this action are analyzed. 

The parameters of the RBC/UKQCD gauge field configurations used here are given in Table llj All ensembles have 
L5 = 16. There are ensembles with two different values of the bare gauge coupling, here given as /? = 6/5^. The 
two values /? = 2.13 and f3 = 2.25 correspond to lattice spacings of a w 0.11 fm and a ~ 0.08 fm, respectively (see 



Sec. IIIB2). The box sizes at the coarser lattice spacing are about 1.8 fm {L ~ 16) and 2.7 fm (L — 24); for the finer 
lattice spacing the box size is about 2.7 fm (L — 32). 

The bottom quark is implemented with lattice NRQCD (SUTU]. The Euclidean action has the form 
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phys.) 


16^ X 32 


2.13 


0.01 


0.04 


2.536 


0.8439 


500- 


4010, 10 


0.1117(33) 


0.436(14) 


2.469(72) 


16^ X 32 


2.13 


0.02 


0.04 


2.536 


0.8433 


500- 


4040, 10 


0.1170(32) 


0.548(16) 


2.604(75) 


16^ X 32 


2.13 


0.03 


0.04 


2.536 


0.8428 


500- 


7600, 10 


0.1195(24) 


0.639(14) 


2.689(56) 


24^ X 64 


2.13 


0.005 


0.04 


2.3, 2.536, 2.7 


0.8439 


915- 


8665, 25 


0.1119(17) 


0.3377(54) 


2.487(39) 


24^ X 64 


2.13 


0.01 


0.04 


2.3, 2.536, 2.7 


0.8439 


1475 


- 8525, 25 


0.1139(19) 


0.4194(70) 


2.522(42) 


24^ X 64 


2.13 


0.02 


0.04 


2.3, 2.536, 2.7 


0.8433 


1800 


- 3600, 25 


0.1177(29) 


0.541(14) 


2.622(70) 


24^ X 64 


2.13 


0.03 


0.04 


2.3, 2.536, 2.7 


0.8428 


1275 


- 3050, 25 


0.1196(29) 


0.641(15) 


2.691(66) 


32^ X 64 


2.25 


0.004 


0.03 


1.75, 1.87, 2.05 


0.8609 


580- 


6840, 20 


0.0849(12) 


0.2950(40) 


1.831(25) 


32^ X 64 


2.25 


0.006 


0.03 


1.75, 1.87, 2.05 


0.8608 


552- 


7632, 24 


0.0848(17) 


0.3529(69) 


1.829(36) 


32^ X 64 


2.25 


0.008 


0.03 


1.75, 1.87, 2.05 


0.8608 


540- 


5920, 20 


0.0864(12) 


0.3950(55) 


1.864(27) 



TABLE I: SuiTimary of lattice parameters for the RBC/UKQCD ensembles. The values for the lattice spacing are results of 
this work and are determined from the T(25) — T(1S) energy splitting (see Sec. Ill B 2 1. The bare gauge coupling is given 
as /3 = 6/g^. The pion masses in lattice units were taken from |16l 1241 132| and converted to physical units using the lattice 
spacings given here. The last column gives the value of the bare b quark mass that would yield agreement of the r;f,(lS') kinetic 
mass with experiment (see Appendix [B|. 



where ip is the two-component bottom quark field, and K{t) is given by 
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Here, Hq is the order-u^ term. 



and 5H contains higher-order corrections, 



i?n = 



A(2) 
2mb ' 



(3) 



SH 



'Cl 



(A(2)) 
8ml 



Cl 



J5_ 
8to2 



V E E V 



'Cz 



+C5 



8m2 



Vx£;-£;xV-C4 -^^— a ■ 



B 



2A(4) 



fA(2)) 



24TOfc 

9 (^{2 



C6 



16n rri^ 



c7/^(a(2),^.b| 

8rn? I J 



cs 



64™^ 



{a(2),^.( 



V X E - E xV 



)} 



cg 



8ml 



(T ■ (E X E) 



(4) 



Note that antiquark propagators can be obtained from quark propagators calculated with the action M through 
Hermitian conjugation. 



Above, A(2) = J2t 



=1^. 
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and A(4) = J2t 



-1^, 
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where V/f denotes the p-ih order symmetric and maximally 



local covariant lattice derivative in //-direction. All derivatives are understood to act on all quantities to their right. 



The chromo-electric and chromo-magnetic fields are defined as Ej = Fjq, Bj 



'2^jkl^kl 



Fki . The terms with coefficients 



Cl to C4 in (Wl) are the order-w^ corrections, while the terms with coefficients C7 to cg are the spin-dependent order-w^ 
corrections (note that spin- independent order-w^ terms are not included). The terms with coefficients C5 and cg are 
spatial and temporal discretization corrections for Hq. Quantities with a tilde also include discretization corrections: 
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In (Im, V)j ''' is a second-order adjoint derivative (which acts only on F^jy) and F^i, is the standard cloverleaf lattice 
gluon field strength. 

At tree- level, the coefficients c^ are equal to 1. The action is tadpole-improved [33] using the mean link in Landau 
gauge, uql • Using the mean link instead of the plaquette for tadpole improvement leads to a better scaling behavior 

m- 

In this work, calculations were performed either with cy = cg = cg = or with cj — cg — cg = 1. These two actions 
will be referred to as the w* action and the v^ action, respectively. The v"^ action is identical to the action used in [TB] . 
The stability parameter n in (pi) was always set to n = 2. Calculations were performed for multiple values of the bare 
6-quark mass a m;,, as shown in Tablell] The spin-dependent energy splittings, which show significant a?7it,-dependence 



(see Appendix C5| were then interpolated to amj^ ^^ , where am\^ ^^'' is the value of the bare 6-quark mass that 



would yield agreement of the ?7t,(15) kinetic mass with experiment (see Appendix [Bj). 

B. Calculation and fitting of tw^o-point functions 

The interpolating fields for the bottomonium two-point functions used here are the same as in |23| , except that the 
cut-off radius was chosen differently (equal to L/2). For a ~ 0.08 fm the smearing parameters in lattice units were 
rescaled from those used at a « 0.11 fm so that the smearing functions in physical units remain the same. 

The methods used here for fitting the two-point functions are also the same as in |23| , i.e. multi-exponential Bayesian 
matrix-correlator fitting combined with statistical bootstrap (suitably modified for Bayesian fitting) [3S]. As in j23| . 
bottomonium two-point functions were calculated for 32 different source locations on each gauge field configuration to 
increase statistics. Note that in ,23j the data were averaged over those source locations prior to the analysis. However, 
the reduced sample size can lead to overestimates of errors due to poorly determined data correlation matrices. As 
shown in Appendix [A] for the i = 24 and L = 32 ensembles, the bottomonium data from the 32 source locations are 
in fact sufficiently independent. Therefore, in the present work the data correlation matrices are calculated with the 
unblocked data sets for the L ~ 2A and L — 32 ensembles. For the L = 16 ensembles, some autocorrelations between 
the data from different source locations were seen, and therefore binning over source locations was performed. 

C. Chiral extrapolations 

The dependence of the bottomonium energy splittings on the light-quark masses is expected to be weak. Therefore, 
chiral extrapolations to the physical pion mass are performed linearly in mj. Before chiral extrapolation, the energy 
splittings are converted to physical units using the lattice spacing determinations on the individual ensembles. In 
addition, the spin-splittings are interpolated to the physical 6-quark mass on each individual ensemble (see Appendices 



C5jand|Bj). 

For a given energy splitting _E(7tt,^, a) that depends on the pion mass m,^ and the lattice spacing a, the chiral 
extrapolation of the data from the L = 32 ensembles with lattice spacing ai w 0.08 fm and the L = 24 ensembles 
with lattice spacing a2 ~ 0.11 fm is performed simultaneously, using the functional form 



E{m^, ai) = E{0,ai) + A', 

E{m^,a2) = E{{),a2)+Ami. (7) 



2 



The three unconstrained fit parameters are £^(0,ai), £(0,02) and A. This allows for an arbitrary dependence of the 
energy splitting E{m.,^, a) on the lattice spacing a. Higher-order terms depending on both a and m^ are neglected. 
As will be shown later, the a-dependence in most quantities is very weak, and therefore these higher-order effects are 
small. No extrapolation in a is performed here, since one can not take the continuum limit with NRQCD. Instead, 
discretization errors can be estimated from the difference of the two results -B(0, ai) and E{0, a2). The simultaneous 
chiral extrapolation was found to significantly improve the statistical accuracy of E(0, oi) from the L = 32 ensembles, 
due to the wider range in m'^ on the L — 2A ensembles. 

The data from the L = 16 ensembles, which have a smaller volume, are extrapolated independently, to allow for an 
arbitrary volume-dependence of the coefficients A (the L = 24 and i = 32 ensembles have approximately the same 
spatial volume in physical units). 

When performing the extrapolations, the statistical uncertainties in the physical pion mass values on the individual 
ensembles (due to the scale uncertainty) are taken into account, by making the pion masses themselves also parameters 
of the fit, constrained with Gaussian priors. The central values and widths of these priors were set equal to the values 
and uncertainties of the pion masses given in Table [l] 



III. RESULTS 

A. Speed of light 

In the continuum, the energy of a particle with mass M and three-momentum p satisfies the relationship E^ = 
M^ +P^ (in the units used here). Equivalently, one has 

^^ ^ L (8) 

To study deviations from the relativistic continuum energy-momentum relationship on the lattice, in the following 
the square of the "speed of light" ([8| will be calculated for the -qbi^S) meson. Note that due to the use of NRQCD, 
energies extracted from fits to bottomonium two-point functions are shifted by approximately —2mb. This shift does 
not affect energy differences. To obtain the full mass of a bottomonium state, one can calculate the kinetic mass, 
defined as 

"'" 2 [E{p) - E{Q)] ^ ' 

Equation (|9| will be equal to the physical mass if the relativistic energy-momentum relationship is satisfied up to 
a constant shift; with lattice NRQCD this is not exact and Mkin will depend slightly on the momentum p. On a 
lattice with L points in the spatial directions and periodic boundary conditions, the momentum can have the values 
p = n ■ 27r/(aL) where n = (ni,n2,'^3) with rii £ Z and —L/2 < rii < L/2. We therefore define the square of the 
speed of light as 

2 [E{p) - E{0) + A4i„i]^ - M^,^^, 
c = -^ , (10) 

where Mkin,i denotes the kinetic mass calculated with 'n? = 1. 

Deviations of c^ from 1 can be caused by missing relativistic and radiative corrections in the NRQCD action (mainly 
in the coefficients Ci, C5 and cg) and by gluon discretization errors. 

The numerical results for c^, calculated using the v'^ action from the rib{lS) energies, are given in Table O and 
plotted in Fig. fl] The results at a « 0.11 fm given here have smaller statistical errors than the previous ones in j23| . 
At a « 0.11 fm, a very small deviation of c^ from 1, at the level of about 0.1% (1.5 a), can now be resolved for n^ = 2 
and n? = 3. At a « 0.08 fm, this deviation goes away, indicating that discretization errors are indeed smaller at the 
finer lattice spacing. 

The statistical error in c^ grows with n^. Notice however that even at n? = 12, which corresponds to a meson 
momentum of about 1.6 GeV, the deviation of c^ from 1 is found to be less than 0.3% at a w 0.11 fm and less than 
0.4% at a « 0.08 fm. This demonstrates that the lattice NRQCD in combination with the Iwasaki gluon action works 
very well for bottomonium. 



n^ 


c^ (asiO.ll fm) 


c^ (a « 0.08 fm) 
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1.00070(47) 


1.00012(55) 


3 


1.00134(85) 


1.00022(92) 
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0.9987(12) 


0.9987(15) 


5 


0.9998(15) 


0.9991(17) 


6 


1.0008(18) 


0.9995(20) 


8 


1.0005(22) 


0.9995(28) 


9 


1.0011(25) 


1.0001(31) 


12 


1.0015(33) 


1.0009(44) 



TABLE II: Square of the speed of light, calculated with the w^ action. For n^ — 9, the components are n = (2, 2, 1) 
and octahedral transformations thereof. The data shown are from the L = 24 ensemble (a ~ 0.11 fm) with ami ~ 0.005, 
am.b = 2.536 and from the L = 32 ensemble (a ~ 0.08 fm) with am,i = 0.004, ami, = 1.87. 
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FIG. 1: Square of the speed of light, calculated with the u* action (points from the two different lattice spacings are offset 
horizontally for legibility). The points at n^ = 12 correspond to a meson momentum of about 1.6 GeV. Note that the very 
small deviation of c^ from 1 at a « 0.11 fm, seen for n^ = 2 and n^ = 3, disappears at a ~ 0.08 fm. 



B. Radial/orbital splittings and the lattice spacing 



Estimates of systematic errors 



The radial and orbital bottomonium energy splittings are calculated in this work with the v'^ action on all ensembles. 
Radial and orbital energy splittings are an effect of order v'^, and therefore the relative error in the radial and orbital 
energy splittings due to the missing v^ corrections is of order v ~ 0.01. Further systematic errors are introduced by 
the missing radiative corrections to the v'^ terms, which leads to a relative error of order Ugv'^ ~ 0.02 (here and in the 
foUowing, Us « 0.2 is used). 

The dominant discretization errors in radial and orbital energy splittings are expected to be caused by missing 
radiative corrections to the couplings C5 and Cg in (El), and by discretization errors in the gluon action. 

The relative errors caused by the radiative corrections to C5 and cg can be estimated as follows: by replacing every 
derivative with a factor of mi,v, we see that the terms with C5 and cg in (H) should be of order a^TO^w'*/(24?7if,) and 
a to|u^/(16 n TO^), respectively. The radiative corrections should be of order as times these estimates. The relative 
error in radial and orbital energy splittings is obtained through dividing by the estimate TO^w^/(2m(,) of Hq. The 
relative error due to the missing radiative corrections to C5 is then 



2 ,,2 
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which is about 1% at a •• 
to Cg is 



0.11 fm and 0.6% at a sa 0.08 fm. The relative error due to the missing radiative corrections 



as a rUb v 
8n 



(12) 



which is about 0.3% at a « 0.11 fm and 0.2% at a sa 0.08 fm. More sophisticated estimates of these errors, using 
wave functions from a potential model, have been made in |18| . 

The discretization errors caused by the gluon action are discussed in detail in Appendix ID] Estimates using tree- 



level perturbation theory are derived in Sec. D 1 The values for the Iwasaki action at the lattice spacings used here 



can be found in Table XXX The nonperturbative results presented in Sees. IIIB2 and D2 indicate that the errors 
may actually be smaller than the tree-level estimates. 



2. Results 



The results for the radial and orbital energy splittings in lattice units at ami, = 2.536 (for a « 0.11 fm) and 
arub = 1.87 (for a w 0.08 fm) are given in Appendix C 1 The spin-averaged masses are denoted with a bar, and are 



defined as 



Ej(2J+1) ■ 



(13) 



In Appendix |C 2[ results for the radial and orbital energy splittings for multiple values of arrib, varying by about 
15%, are given. As can be seen there, the dependence on anib is very weak. The change when interpolating from 
ami, = 2.536 at a « 0.11 fm or ami, — 1.87 at a sa 0.08 fm to the physical values of ami, given in Table IT] would be 
much smaller than the statistical errors. 

The inverse lattice spacings of the gauge field ensembles are determined from the T(2S') — T(15) splitting, dividing 
the experimental value from [36] by the dimensionless lattice value. The 25*— IS* spl itting is expected to have smaller 
systematic errors than the IP — IS splitting [TB]. In particular, as shown in Sec. ]D lb[ the gluonic discretization 
errors partially cancel in the 25" — 15 splitting. For the Iwasaki action and at tree-level, the remaining gluon errors 
in the 2S — IS splitting are estimated to be about 2.6% at a « 0.11 fm and 1.6% at a « 0.08 fm. 

Results for the lattice spacings of all ensembles are given in Tables |III[ |IV[ an d \V\ The errors shown there are 
statistical/fitting only. For comparison, results from both the T(2S') — T(IS') and l-^P — T(IS') splitting are given. 
They are found to be mostly consistent within the statistical errors here. Note that in the quenched approximation, 
the ratio of the 25* — IS* and IP — IS* splittings was previously found to be in disagreement with experiment ^18J . 

The lattice spacings obtained here are seen to become slightly finer as the sea quark mass is reduced (in |23j . 
this dependence was hidden by the larger statistical errors). This behavior is in fact expected here, since in the 
RBC/UKQCD ensembles, the bare gauge coupling is kept constant when varying the sea-quark masses (see Table ]l|. 
In contrast to this, the MILC collaboration decreases /3 ex 1/g^ when decreasing the sea-quark masses [13], such that 
the lattice spacing remains approximately constant |18) . 
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0.01 



anil 



0.02 
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0.03 






1.766(52) 
1.718(16) 



1.687(46) 
1.678(13) 



1.651(33) 
1.661(10) 



am; =0.005 am; =0.01 am; =0.02 am; =0.03 






-IS 



1.763(27) 1.732(28) 1.676(42) 1.650(39) 
1.742(14) 1.703(12) 1.680(27) 1.670(33) 



TABLE III: Inverse lattice spacings of the L 
in GeV. 



16 ensembles TABLE IV: Inverse lattice spacings of the L • 

in GeV. 



24 ensembles 



ami = 0.004 ami = 0.006 ami = 0.008 



^2S- 



IS 
IS 



2.325(32) 
2.305(24) 



2.328(45) 
2.329(23) 



2.285(32) 
2.328(23) 



TABLE V: Inverse lattice spacings of the L = 32 ensembles in GeV. 



The lattice spacings from the 25—15' splittings on the individual ensembles were then used to convert the results 
for the other radial and orbital energy splittings to physical units. The conversion was performed using the bootstrap 
method to take into account correlations between the 25 — 15 splitting and the other splittings. 

Finally, the results in physical units were extrapolated to the physical pion mass, as shown in Fig. [2] The data 
from the L = 24 and L = 32 ensembles (both have a box size of about 2.7 fm) were extrapolated simultaneously using 
the function Q. The data from the L — 16 ensembles were treated independently, since the dependence on the pion 
mass may be different for the smaller box size of about 1.8 fm. 
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FIG. 2: Chiral extrapolation of radial and orbital energy splittings (calculated with the d* action). Extrapolated points are 
offset horizontally for legibility. 



T(35') - T(IS') 
]AP-T(1S') 

Fp-Is 
Wp-r{is) 

WP-TS 

T2(1D)-T(15') 



L = 16, a ^ 0.11 fm 
0.98(13) 
0.462(24) 
0.476(24) 
0.414(61) 
0.865(70) 
0.878(70) 
0.717(39) 



L = 24, a ^ 0.11 fm 
0.883(52) 
0.4524(78) 
0.4657(79) 
0.351(21) 
0.802(23) 
0.816(23) 
0.709(14) 



L = 32, a ^ 0.08 fm 
0.900(34) 
0.4428(57) 
0.4571(59) 
0.362(15) 
0.803(17) 
0.817(17) 
0.7060(98) 



Experiment 
0.89490(56) 
0.43957(37) 
0.45690(79) 
0.36033(45) 
0.79990(47) 
0.81722(83) 
0.70340(87) 



TABLE VI: Chirally extrapolated results for the radial and orbital energy splittings in GeV, calculated with the v'^ action. 
The errors on the lattice results shown here are statistical/fitting/scale setting only; for a discussion of systematic errors, see 
Sec. ~ 



IIIBl 



The experimental values are from [TH3] (for the 'rit{lS), which enters in IS), |37l I38j (for the T2(1D)) and from 
the Particle Data Group i36j (for all other states) 
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FIG. 3: Chirally extrapolated results for the radial and orbital energy splittings, calculated with the u* action. The T(15') and 
T(2S) masses are used as input in the lattice calculation. The errors on the lattice results shown here are statistical/fitting/scale 
setting only. 



The numerical results at the physical pion mass are given in Table VI In addition, the energies of the radial and 
orbital excitations at the physical pion mass are plotted in Fig. [3| The extrapolated data from the L = 16 ensembles 
have significantly larger statistical errors than the other data. This is expected for the following two main reasons: 
first, the four-dimensional volume of these ensembles is about 7 times smaller than that of the L — 24 ensembles, 
providing less information. Second, the lowest pion mass available at L = 16 is larger than on the other two ensembles, 
requiring more extrapolation. Within the statistical errors, no finite-volume effects are seen (also at the individual 
values of the quark masses, where the results are more precise; see Fig. p|. Finite- volume effects in bottomonium have 
been studied in detail using a lattice potential model in [3^. At a box size of 1.8 fm (corresponding to the L = 16 
lattices here) the most significant finite-size effects were found in the 35* energy, which was shifted by about —40 MeV 
compared to the infinite-volume energy. At a size of 2.7 fm, this shift was found to be negligible. 

The results for the radial and orbital energy splittings obtained here from the L — 24 and L = 32 ensemble show 
very little dependence on the lattice spacing. When going from a « 0.11 fm to a « 0.08 fm, the IP — IS* splitting 
changes by about 2%. However, this change is only 1 standard deviation and could also be caused by a statistical 
fiuctuation. All results at a w 0.08 fm, where discretization errors are expected to be the smallest, are in excellent 
agreement with experiment. The 1-^P — IS" and T2(l-D) — T(IS') splittings at a w 0.08 fm have statistical errors 
of only about 1.3%, and they agree fully with experiment. This indicates that the systematic errors are very small; 
smaller than the tree-level estimates of gluon-discretization errors given in Table IXXX] would suggest. 
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C. Spin-dependent splittings 



1. Discussion of systematic errors 
Spin-dependent energy splittings first arise at order w^ in the nonrelativistic expansion, via the terms with coefficients 



C3 and C4 in Eq. Hi. For the order-w action, it is expected that at lowest order the S-wa,ve hyperfine splittings 



and the P-wave tensor splittings 



T{nS) - rib{nS) 



are proportional to c1 and independent of C3 , while the P-wave spin-orbit splittings 



^XboinP) + 3X61 (nP) - Xb2{nP) 

C3, while the P-wave spin-orbit sp 

2xboinP) - 3x61 (nP) + 5X62 (nP) 



(14) 

(15) 

(16) 



are proportional to C3 and independent of C4 |18j . In this work, the coefficients Ci are set to their tree-level values, 
Ci = 1. Therefore, spin-dependent energy splittings calculated directly will have systematic errors of order a^ ~ 0.2. 
However, these unknown radiative corrections are expected to cancel in ratios of quantities with equal dependence 
on the couplings Ci . Nonperturbati ve re sults for the dependence of various spin splittings and ratios of spin splittings 
on C3 and C4 are given in Appendix C 4 for both the v'^ and v^ actions. As can be seen there, ratios of hyperfine and 
tensor splittings are indeed independent of both C3 and C4 to a very good approximation. 

Spin-splittings calculated with the u^ action will also have relativistic errors of order u^ sa 0.1 due to the missing v^ 
corrections. Therefore, the spin-dependent t;^-corrections, given by the terms with coefficients cy, cg and cg in Eq. (|4|, 
are included in this work. Relativistic corrections for spin splittings calculated with the v^ action are then expected 
to be of order v'^ « 0.01 due to the missing order-w^ terms. Missing radiative corrections to the order-w^ terms lead 
to additional systematic errors of order agv"^ « 0.02. 

The terms with coefficients C3 and C4 in Eq. Q include the tree-level discretization corrections via (l5| and ^: 



2 



-y^jMc^i I Vfe 



7(2, ad) 



F, 



10 



7(2, ad) 



Fu 



„(2,ad) p yji^M) p 



(17) 



a B = 



a. 

2 - 



^jkl Fi 



kl 



12 



CT7 £jkl 



VI^'^'^^P,; 



vp'^^V,, 



However, radiative corrections to the order-a terms in (17 1 and (18 1 are missing. Their size can be estimated using 



(18) 



the NRQCD power-counting rules. We replace every spatial derivative by mbV and every temporal derivative by 
^rribV^ (by the leading-order equations of motion, a temporal derivative is of the order of the nonrelativistic kinetic 
energy [in])- The radiative discretization corrections are of order a^ times the a^ terms. For the spin-orbit splitting 
controlled by (17), the terms with the spatial derivatives in the a^ terms are dominant, and the relative discretization 
errors become 



1 



2 2 2 

asa m^v . 



(19) 



This is about 0.04 at a ft 
and tensor splittings are 



0.11 fm and about 0.02 at a « 0.08 fm. The relative discretization errors in the hyperfine 



222 



(20) 



Here, an additional factor of 2 was introduced to take into account the quadratic dependence of the hyperfine and 
tensor splittings on (18). Equation (20) is equal to about 0.09 at a « 0.11 fm and about 0.05 at a « 0.08 fm. 

In the spin-dependent order-f^ corrections (the terms with coefficients C7, cg and Cg in the action), the lattice 
Laplacian used here does not include discretization corrections. The tree-level corrected Laplacian would be A^^^ = 
A*^^) — (a^/12)A('*). The relative error in the hyperfine and tensor splittings caused by the missing of this correction 
is then of order o?m\v^ /%^ which is only 0.005 at a « 0.08 fm. 

Additional discretization errors in spin splittings may arise from the gluon action. These errors are discussed in 



detail in Sec. D 2 For the S'-wave hyperfine splittings, the relative error caused by discretization errors in the Iwasaki 
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action is estimated to be about 5% at a sa 0.08 fm; for the other spin sphttings no significant ghion discretization 
errors are found. 

Finally, note that the lattice NRQCD calculation performed here does not include the effects of annihilation of the 
b and b. This mainly affects the pseudoscalar states rii,{nS), which can annihilate into two giuons (for the T(nS), at 
least three giuons are required). The annihilation contribution to the rj^ mass can be related to the two-gluon decay 
width r[T]b -^ gg] as follows [4CT1I4T] : 

-J^annihii. [Vb] « ^"^^^ ^ ^ r[% ^ gg]. (21) 

For the rii,{lS) and r]b{2S), the widths were calculated in 42J to be 7 MeV and 3.5 MeV, respectively. This gives 
SEannihii.[ilb{^S)] = —0.7 MeV and (5_Eaimihii. [?7b(2S')] = —0.34 MeV, which is only about 1% of the hyperfine sphttings. 

2. Results 
The results for the spin-dependent energy splittings in lattice units at ami, — 2.536 (for a « 0.11 f m, L = 24) and 



anib = 1.87 (for a ~ 0.08 fm, L = 32) for both the v^ and v^ NRQCD actions are given in Appendix C3 Here, all 
couplings Ci in the action were set to their tree-level values of 1. 

In Appendix C 4 the dependence of the spin splittings on the couplings C3 and C4 in the action Q is studied (on the 
i = 24 ensemble with ami = 0.005). Results are shown both for the w^ and v^ actions. The n aive ex pectations for the 
C3- and C4-dependence of the order-w^ spin splittings were discussed at the beginning of Sec. Ill C 1 The dependence 



of the IP spin-orbit splitting on the coupling C3 appears to be slightly weaker than expected: it changes by only 13% 
(for the v'^ action) or 15% (for the v^ action), when C3 is changed by 20%. Contrary to the naive expectation, the 
IP spin-orbit splitting also shows some dependence on C4: about 6% (for the v'^ action) or 8% (for the v^ action), 
when C4 is varied by 20%. On the other hand, the IP tensor splitting behaves as expected: it shows no significant 
C3-dependence, and the dependence on C4 is consistent with proportionality to c|, both for the w* and v^ actions. The 
results for the 15 hyperfine splitting are also close to these expectations. However, in the IS hyperfine splitting the 
deviations from the naive expectations, while not large in absolute terms, are statistically significant due to the very 
small statistical errors. The most important result from Appendix |C4| is that the ratio of the 25* and IS* hyperfine 
splitting as well as the ratios of the S-wetve hyperfine and the IP tensor splitting show no significant dependence on 
either C4 or C3, and this is true for both the w* and v^ actions. 

Next, in Appendix |C 5[ results for the spin-dependent energy splittings for multiple values of ami,, varying by 
about 15%, are given, and visualized in Fig. [o] As can be seen there, the results for most splittings are compatible 
with a l/rrib-dependence, with the notable exception of the IS* hyperfine splitting. The l/mb-dependence of the spin 
splittings can be understood from the power-counting rules as follo ws: ra dial and orbital energy splittings, which are 



of order mbv^ , are nearly independent of mi,, as shown in Appendix C 2 This implies that 

w^ ex 1/to6. (22) 

Spin-dependent energy splittings are a factor of v^ smaller than radial and orbital energy splittings. Since the latter 
are nearly constant, spin-dependent energy splittings are expected to be proportional to v^ , and hence l/m,b- 

The results for the 15* hyperfine splitting in lattice units have very small statistical errors, and are clearly incom- 
patible with a dependence proportional to l/{am.i,)- However, fits of the form A/{am,i,) + B with a constant term B 
describe the date very well in the range considered here. The fit results A and B for the IS hy perfine splittings on all 



L = 24 and all L = 32 ensembles, for both the w* action and the v^ action are given in Tables XXIII |XXV XXVII 
andExIXl 

To obtain physical results, all spin splittings were then interpolated to the physical &-quark masses given in Appendix 
[B| assuming a l/(amf,) dependence everywhere except for the IS" hyperfine splittings and the ratios. For the \S 
hyperfine splittings and the ratios involving them, the fit results A and B were used in the interpolation. 

The interpolated spin splittings were then converted to physical units using the lattice spacings from the 25* — 15 
splittings on the individual ensembles as obtained in Sec. |IIIB 2| Note that the uncertainty in the lattice spacing 
enters with a factor of 2 here, due to the resulting uncertainty in the bare heavy quark mass and the approximate 
l/{am,b) behavior of spin splittings. 

Finally, simultaneous chiral extrapolations of the data at a w 0.11 fm and a k, 0.08 fm to the physical pion mass, 
using the functional form ([7|), were performed (the data from the v^ and v^ actions were treated independently). 
These chiral extrapolations are visualized in Figs. [4] and [5] 

The numerical results for the spin-splittings at the physical pion mass are given in Table [VIl| and plotted in Fig. |6] 
It can be seen that the results obtained with the v^ and with the v^ action differ significantly. At a w 0.08 fm, the 15 
hyperfine and spin-orbit splitting are reduced by about 20% by the v^ terms, while the IP tensor splitting is reduced 
by 10%. These changes are in line with the estimate of v"^ ~ 0.1. 
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FIG. 4: Chiral extrapolation of spin-dependent splittings from the L = 24 and L = 32 ensembles, part I. Extrapolated points 
are offset horizontally for legibility. 
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FIG. 5: Chiral extrapolation of spin-dependent splittings from the L = 24 and L = 32 ensembles, part II. Extrapolated points 
are offset horizontally for legibility. 
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V action 
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-2x6o(lP) + 3x6i(lP)-Xw(lP) 

ri2S)-r^b{2S) 
-2x6o(lP)+3xw(lP)-Xw(lP) 



55.7(2.0) 


59.0(1.4) 


46.9(1.7) 


48.5(1.1) 


69.3(2.9) 


22.1(4.1) 


27.1(2.6) 


15.8(3.6) 


19.7(2.4) 


- 


20.9(2.3) 


20.9(1.9) 


17.0(2.0) 


16.9(1.7) 


19.43(57) 


28.9(2.1) 


31.3(1.7) 


26.0(1.8) 


27.5(1.4) 


33.34(66) 


0.9(1.4) 


0.3(1.0) 


0.6(1.2) 


0.04(93) 


- 


163(13) 


167(11) 


137(12) 


139.4(9.8) 


163.8(2.6) 


36.6(4.4) 


41.6(3.3) 


34.7(3.9) 


38.1(3.0) 


47.3(1.6) 



0.391(75) 



1.54(18) 



0.57(14) 



0.456(47) 



1.44(13) 



0.622(98) 



0.333(78) 



0.403(52) 



1.36(16) 1.28(12) 

0.43(13) 0.497(87) 



1.467(79) 



TABLE VII: Spin-dependent energy splittings from the L = 24 and L = 32 ensembles, interpolated to the physical b quark 
mass and extrapolated to the physical pion mass. All results in MeV, except for the dimensionless ratios. For the lattice data, 



the errors shown here are statistical/fitting/scale setting only; see Sec. IIIC 1 for a discussion of systematic errors and Table 
VIII for the final results that include estimates of the systematic errors. The experimental value for T(IS') — rjt{lS) is the 
weighted average of the results from [l], [2], and [3]. All other experimental values are from the Particle Data Group |36| . 
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FIG. 6: Chirally extrapolated results for the spin-dependent energy splittings, from the L = 24 and L = 32 ensembles. Left 
panel: IP spin splittings (energies relative to l^P). Right panel: S'-wave hyperfine splittings (energies relative to the T(IS') 
and T(2S) states, respectively). The errors shown are statistical/fitting/scale setting only; see Sec. IIIC l| for a discussion of 
systematic errors and Table [VIIl] for the final results that include estimates of the systematic errors. 

The dependence on the lattice spacing varies between the different quantities. For example, the IS* hyperfine 
splitting calculated with the v^ action increases by 3% when going from a ~ 0.11 fm to a w 0.08 fm, but this change 
is only 0.8 standard deviations. With the u** action, the lattice spacing dependence in the 15 hyperfine splitting 
appears to be stronger, about 6% or 1.4 standard deviations. This size of the a-dependence is in good agreement 
with the estimates of discretization errors obtained in Sec. |IIIC 1[ Note that the dependence on the lattice spacing in 
spin splittings calculated directly (as opposed to the ratios) may be caused both by discretization errors and by the 
amb-dependence of the missing radiative corrections to C3 and C4. When going from a ~ 0.11 fm to a « 0.08 fm, the 
IP tensor splitting appears to change by about 10%, but the effect is less than 1 standard deviation. No a-dependence 
is seen in the IP spin-orbit splitting. Recall that the estimates of discretization errors obtained in Sec. |IIIC 1| are 
indeed smaller by a factor of 2 in the spin-orbit splitting compared to the tensor and hyperfine splittings. In addition, 
in potential models the spin-orbit splitting is not as sensitive to short distances as the hyperfine splitting. 
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The ratios of spin splittings calculated here also show no significant dependence on the lattice spacing: about 0.7 
standard deviations in the ratio of the 2S and IS* hyperfine splittings and about 0.5 standard deviations in the ratios 
of the S-wscve hyperfine and IP tensor splittings. 

The most reliable results for the spin-dependent energy splittings obtained in this work, calculated with the v^ 



action at a « 0.08 fm, are summarized in Table VIII Here, estimates of the systematic errors based on the discussions 
in Scc. |IIICT] and Scc. |D2| are given. The systematic errors in the ratios of the hyperfine and tensor splittings, where 
the unknown radiative corrections to C3 and C4 cancel, are dominated by discretization errors. The systematic error 
in the l^P — hi,{lP) hyperfine splitting is dominated by the unknown radiative corrections of order as. However, the 
absolute systematic error in l^P — hb{lP) is only 0.2 MeV (assuming as ~ 0.2), because the splitting is found to be 
zero within the statistical error of about 1 MeV. The l^P — hi,(lP) splitting vanishes in potential models because the 
wave function at the origin is zero for L ^ 0. 

Using the lattice ratios of the S'-wave hyperfine and the IP tensor splitting, and the experimental result for the IP 
tensor splitting I3S1, the IS' and 25" hyperfine splittings can be calculated in MeV. The 25* hyperfine splitting in MeV 
can be calculated alternatively from the lattice ratio of the 25" and IS hyperfine splitting, and the experimental value 



of the IS* hyperfine splitting [THl]. The results obtained with both methods are shown in Tabic VIII 



This work Experiment 



T(25') - nt{2S) 

r{is)~vb{is) 

r{is)-vb{is) 

-2xboilP) + SxbiilP) - Xb2{lP) 
r{2S)-Vb{2S) 



0.403(52)(27) 

1.28(12)(10) 1.467(79) 

0.497(87) (44) 



-2xbo{lP) + 3xbi(lP) - Xw(lP) 

T(15') - VbilS) 60.3(5.5)(5.0)(2.1) MeV " 69.3(2.9) MeV 



23.5(4.1)(2.1)(0.8) MeV " 
28.0(3.6)(1.9)(1.2) MeV ' 



T(25')- 7,6(25) 

FP - ht{lP) 0.04(93)(20) MeV 



"Using IP tensor splitting from experiment 
Using T(1S) — rib{lS) splitting from experiment 

TABLE VIII: Final results for spin splittings, calculated with the v^ action at a ~ 0.08 fm (L = 32). For the lattice data, the 
first error is statistical/fitting, the second error is an estimate of systematic uncertainties, and the third error (where given) is 
experimental. The experimentai vafue for the IS hyperfine splitting is the weighted average of the results from [T], [2], and [3]; 
the experimental value for the IP tensor splitting is calculated using the I'^P masses from the Particle Data Group [36) . 

Using the T(IS'), T{2S) and I'^P masses from experiment [36,, the absolute masses of the r]i,{lS), riii{2S) and 
/if, (IP) mesons can then be calculated. This gives 

M[tji,{1S)] = 9.4000(55)(50)(21) GeV, 

M[Tjbi2S)]'' = 9.9998(41)(21)(9) GeV, 

M[r]b(2S)]'' = 9.9953(36)(19)(12) GeV, 

M[hbilS)] = 9.89983(93)(20)(27) GeV, (23) 

where the first error is statistical/fitting, the second error is systematic, and the third error is experimental. For the 
i]i,{2S) mass, the results from both methods as discussed above are given. 
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IV. CONCLUSIONS 



In this paper, a high-precision calculation of the bottomonium spectrum in lattice QCD with 2 + 1 flavors of 
dynamical light quarks was presented. One important improvement over j23j was the inclusion of a finer lattice 
spacing, giving better control of discretization errors. The dependence of the results on the lattice spacing was seen 
to be weak, and consistent with the estimates of NRQCD discretization errors based on power counting. At a « 0.08 
fm, the radial and orbital energy splittings were found to be in excellent agreement with experiment, within statistical 
errors as small as 1.3% (see Fig.|3|. In addition, the square of the "speed of light", a quantity used on the lattice to 
measure deviations from the relativistic continuum energy-momentum relationship, was found to be compatible with 
1 within statistical errors smaller than 0.4% for bottomonium momenta up to 1.6 GeV (see Fig. fTl). These results 
provide valuable tests of the lattice methods used here: NRQCD for the b quarks, the Iwasaki action for the gluons, 
and the domain wall action for the sea-quarks. The discretization errors associated with the gluon action were studied 
further using a lattice potential model based on tree-level perturbation theory, and by repeating the nonperturbative 
calculations on MILC gauge field ensembles generated with the Liischer-Weisz action. These tests show that the 
Iwasaki action works well at the lattice spacings considered here. 

The focus of this work was the accurate calculation of spin splittings. To this end, ratios of hyperfine and tensor 
splittings were calculated, in which the unknown radiative corrections to the leading spin-dependent terms in the 
NRQCD action cancel. This cancellation was confirmed here directly through numerical calculation of these ratios 
with different values of the spin-dependent couplings in the action. Furthermore, systematic errors from relativistic 
corrections were reduced by the inclusion of the spin-dependent order-u^ terms in the NRQCD action. The results in 
Table [VIII| are considerably more precise than those from previous lattice computations. For example, the ratio of the 
25' and IS* hyperfine splittings is predicted here to be 0.403±0.052 stat± 0.027 syst (the result from 181 is 0.5±0.3 stat)- 



By the criterion of Ref. [43], the results (23) for the ?7&(25') mass obtained here are now the most accurate predictions 
of a gold-plated hadron mass from lattice QCD to date. The prediction of the hi,{lP) mass appears to be even more 
accurate, but note that it is obtained by subtracting from the experimental result for the l^P mass a splitting that 
is zero within the statistical errors. 

The result for the IS hyperfine splitting obtained here is 60.3 ± 5.5 stat i 5.0 gyst i 2.1 exp MeV. This is consistent 
with the value of 54 ± 12 MeV calculated with the Fermilab method in [T3]. It is also in excellent agreement with the 
prediction of 60 MeV obtained in [44] using a relativistic quark model, and with the result of 58 ± 1 MeV from [45] 
for Hf = 3. 

The 15 hyperfine splitting calculated here is only about 1 standard deviation below the weighted average of the 
experimental results from [IH3]. The splitting obtained here is larger than many results from pcrturbative QCD 
[IHS]. Penin argues in [6] using continuum perturbation theory (where in fact lattice perturbation theory should be 
used) that the inclusion of radiative corrections in the NRQCD action could reduce the lattice value for the hyperfine 
splitting by about 20 MeV and bring it in agreement with pcrturbative QCD. This statement does not apply to the 
result obtained here, where the hyperfine splitting is calculated from the ratio to the P-wave tensor splitting so that 
the radiative corrections cancel. Interestingly, it is noted in [46] that the pcrturbative prediction for the 15 hyperfine 
splitting increases significantly and becomes consistent with experiment when the delta function terms in the potential 
are not softened. 
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Appendix A: Analysis of autocorrelations 

In this work, bottomonium two-point functions were calculated for 32 different source locations spread evenly 
across the lattice on each gauge field configuration in order to increase statistics. The question is whether data 
from different source locations are statistically independent, and also whether the data from successive (in molecular 
dynamics time) gauge field configurations are statistically independent. Possible autocorrelations in the data can 
be reduced by binning, that is, by averaging the data within blocks of some size B prior to the further statistical 
analysis. Increases in the statistical errors in an observable under binning of the data would indicate the presence of 
autocorrelations. 

Performing the binning analysis for the energies obtained from the matrix fits used here with their large number 
of degrees of freedom is problematic due to spurious finite-sample-size effects for the estimates of the data correlation 
matrix when the bin size becomes too large. Therefore, in the following analysis, the statistical errors in the two-point 
functions themselves are considered instead. Figure ITJ shows the errors in the diagonal T(IS'), T(2S) and T(3S') 
two-point functions for a given source-sink separation i, versus the bin size B (all errors relative to the corresponding 
error at _B = 1). To the left of the vertical dashed lines in the graphs, binning is performed over neighboring source 
locations only {B — 1, 2,4, ...,32). To the right of the vertical dashed lines, binning is performed over neighboring 
gauge configurations {B = 64, 96, 128, ..., 256). Here, "neighboring" gauge configurations are separated by the step 
sizes given in Tablejl] As can be seen in Fig.JTJ for the L = 16 ensemble, which has a box size of about 1.8 fm, significant 
autocorrelations between the data from the different source locations can be seen at short time separations for the 
interpolating fields optimized for the excited states (T(2S') and T(3S')). The stronger autocorrelations for excited 
states compared to ground states can be explained by the larger physical size of the excited states. No significant 
autocorrelations are seen in molecular dynamics time. Note that the source locations were always shifted randomly 
from configuration to configuration in this work. For the L — 2A and L — 32 ensembles, which have a box size of 
about 2.7 fm, no significant autocorrelations are seen either between source locations or in molecular dynamics time. 
The same qualitative behavior was found for other bottomonium two-point functions. 
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FIG. 7: Analysis of autocorrelations in the two-point functions of the T(15'), T{2S) and T(35') interpolating fields. The data 
are from the L — 16 ensemble with ami = 0.01, the L — 24 ensemble with ami = 0.005, and the L = 32 ensemble with 
ami = 0.004. 
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Appendix B: Tuning of the bare 6 quark mass 



Results for the kinetic mass of the ?7b(15) meson in lattice units, calculated from Eq. dol) with the smallest possible 
magnitude of the lattice momentum |p| = 1 • 27r/(aL), are given in Tables IX and Ix] for the L — 24 and L — 32 
ensembles, respectively. Note that aMkin is nearly independent of p, as demonstrated by the calculation of the "speed 
of light" in Sec. Ill A As can be seen in Fig. |8] the dependence of aMkin on ami, is consistent with the linear function 
aMkin = A ■ ami, + -B in the ranges considered. The fit results for A and B are also shown in the tables. Using 
these fit results, and the lattice spacings from the T(2S') — T(IS') splitting given in Tables IV andlvj the "physical" 

values of the bare b quark mass amf^ ^^'' were then determined such that Mkin agrees with the experimental value of 
9.3910(29) GeV [IHS]. _ 

Note that using the spin-averaged 15' kinetic mass instead of the 77b (15) kinetic mass for the tuning gives values 

of amj^ '^^ that are about 1% larger. The resulting shifts in the spin splittings (in physical units) would be smaller 
than the statistical errors obtained here. For the ratios of spin splittings the shifts would be only about 0.05 standard 
deviations, which is negligible. 



The tuning of the 5-quark mass was performed here only for the v action, and the values of am 



(phys.) 



obtained with 



this action were then used also for the calculations with the v^ action. The v^ action employed here does not include 
the spin-independent order-ii^ terms, and therefore a calculation of Mkin with this action would not be complete to 
this order. Ignoring this issue and doing the tuning for the v^ action was found to increase am,j^ ^^'' by about 2% 
relative to the values obtained for the w^ action using the 15 kinetic mass. Again, the effect of this shift on the ratios 
of spin splittings would be negligible. For the v^ action the values of am,j^ ^^ calculated using the kinetic masses of 
the 77b(15), the spin average 15, and the T(15) were found to be in agreement. 



anib = 2.3 



amb — 2.536 



anib — 2.7 



Fit result 



(phys.) 



ami = 0.005 
ami = 0.01 
ami = 0.02 
ami = 0.03 



4.965(11) 
4.986(13) 
4.979(45) 
4.933(29) 



5.414(22) 
5.447(15) 
5.443(50) 
5.391(31) 



5.743(13) A = 1.944(12), B = 0.494(27) 2.487(39) 

5.768(17) A = 1.958(18), B = 0.483(44) 2.522(42) 

5.763(54) A = 1.967(56), B = 0.46(12) 2.622(70) 

5.712(32) A = 1.947(38), B = 0.454(95) 2.691(66) 



TABLE IX: Heavy-quark mass dependence of the kinetic mass of the 776(15) meson in lattice units, calculated with the u'' 
action on the L = 24 ensembles. Also shown are the results of correlated fits using the functional form aMkin — A ■ amb -\- B, 
and the value of anij, that would yield agreement of the rib{lS) kinetic mass with experiment. 



anib ~ 1.75 



anib = 1.87 



am,b = 2.05 



Fit result 



(phys.) 



am; =0.004 3.882(10) 4.112(11) 4.458(11) A = 1.914(16), B = 0.534(29) 1.831(25) 

ami = 0.006 3.8823(85) 4.1118(93) 4.4560(98) A = 1.912(10), B = 0.536(19) 1.829(36) 

ami = 0.008 3.889(13) 4.122(14) 4.468(15) A = 1.927(22), B = 0.518(40) 1.864(27) 



TABLE X: Heavy-quark mass dependence of the kinetic mass of the r;t,(lS) meson in lattice units, calculated with the ii* action 
on the L = 32 ensembles. Also shown are the results of correlated fits using the functional form aMkin = A ■ arub + B, and the 
value of anib that would yield agreement of the rjb{lS) kinetic mass with experiment. 
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FIG. 8: Kinetic mass of the 7^5 (15*) meson plotted as a function of the bare heavy quark mass (both in lattice units). The lines 
and error bands are from correlated fits using the functional form aMkin ~ A ■ arrib + B. Left panel; L = 24, ami ~ 0.005, right 
panel: L = 32, arm = 0.004. 
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Appendix C: Results in lattice units 
1. Radial and orbital energy splittings 





ami = 0.01 




ami = 0.02 




ami = 0.03 


T(25') - T(1S) 


0.3187(93) 




0.3337(91) 




0.3409(68) 


2S-1S 


0.3233(84) 




0.3372(84) 




0.3448(70) 


T(3S') - T(IS') 


0.536(42) 




0.545(37) 




0.549(26) 


l3p-r{lS) 


0.2559(24) 




0.2620(20) 




0.2646(15) 


lip -IS 


0.2634(24) 




0.2697(20) 




0.2725(15) 


23p^l3p 


0.226(21) 




0.244(16) 




0.236(16) 


2'iP-T{lS) 


0.478(23) 




0.506(17) 




0.501(16) 


23P - IS 


0.485(23) 




0.513(17) 




0.509(16) 


T2(li5) - T(1S) 


0.4052(84) 




0.4254(60) 




0.4318(42) 


TABLE XI: Radial and orbital energy splittings in 


lattice units, 


calculated with the ii* actior 


L on the L = 16 ensembles, for 


anib — 2.536. 














ami = 0.005 


ami = 0.01 




ami =0.02 


ami = 0.03 


T(2S')-T(15) 


0.3193(49) 


0.3250(53) 




0.3359(84) 


0.3413(82) 


2S-1S 


0.3236(46) 


0.3291(51) 




0.3395(87) 


0.3453(81) 


T{3S) - r(is) 


0.497(22) 


0.542(23) 




0.540(53) 


0.573(51) 


13P-T(1S') 


0.2523(20) 


0.2580(19) 




0.2617(42) 


0.2632(51) 


13P-1S 


0.2598(20) 


0.2656(19) 




0.2695(42) 


0.2710(51) 


23p_l3p 


0.207(12) 


0.2240(79) 




0.246(13) 


0.263(19) 


23p-r{is) 


0.460(14) 


0.4820(84) 




0.508(15) 


0.526(21) 


2'iP~lS 


0.467(14) 


0.4896(84) 




0.516(15) 


0.534(21) 


r2{lD)-T{lS) 


0.3998(54) 


0.4140(49) 




0.4330(82) 


0.430(11) 



TABLE XII: Radial and orbital energy splittings in lattice units, calculated with the v action on the L = 24 ensembles, for 
amt — 2.536. 



ami = 0.004 



ami = 0.006 



ami = 0.008 



T(2S')-T(15) 

2S-TS 

r{3S) ~ T(IS') 

1AP-T(15) 

I3P-IS 

WP-WP 

Wp-r{is) 
Wp-Is 

T2(1D) - T(15') 



0.2421(33) 

0.2454(33) 

0.394(15) 

0.1907(20) 

0.1969(20) 

0.1629(94) 

0.3519(94) 

0.3580(94) 

0.3051(40) 



0.2418(47) 

0.2452(45) 

0.395(16) 

0.1888(19) 

0.1949(19) 

0.1649(85) 

0.3533(94) 

0.3594(94) 

0.3045(56) 



0.2464(34) 

0.2493(34) 

0.401(16) 

0.1888(19) 

0.1950(19) 

0.170(11) 

0.359(11) 

0.365(11) 

0.3088(54) 



TABLE XIII: Radial and orbital energy splittings in lattice units, calculated with the «■* action on the L = 32 ensembles, for 
arub — 1.87. 



2. Heavy-quark mass dependence of radial and orbital energy splittings 
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T(2S')-T(15) 

2S-TS 

1AP-T(15) 

Fp-Is 

WP-WP 

Wp-r{is) 
Wp-Is 

T2(1D)-T(15') 



amt — 2.3 

0.3211(50) 

0.3257(47) 

0.2509(21) 

0.2590(21) 

0.213(15) 

0.464(16) 

0.472(16) 

0.4002(52) 



arrib = 2.536 

0.3193(49) 

0.3236(46) 

0.2523(20) 

0.2598(20) 

0.207(12) 

0.460(14) 

0.467(14) 

0.3998(54) 



arrit = 2.7 

0.3184(49) 

0.3225(48) 

0.2536(19) 

0.2608(19) 

0.204(11) 

0.458(12) 

0.465(12) 

0.4001(57) 



TABLE XIV: Heavy-quark mass dependence of radial and orbital energy splittings in lattice units, calculated with the w* action 
on the L = 24 ensemble with ami ~ 0.005. 



anih = 1.75 



arrib ~ 1.87 



arub = 2.05 



T(2S') - T(IS') 
2S-TS 
I3P-T(1S') 
RP-TS 

Wp-T^p 

¥p-r{is) 

Wp-Ts 

r2{lD)-T{lS) 



0.2422(31) 

0.2456(32) 

0.1901(22) 

0.1965(22) 

0.1645(99) 

0.353(10) 

0.359(10) 

0.3048(39) 



0.2421(33) 
0.2454(33) 
0.1907(20) 
0.1969(20) 
0.1629(94) 
0.3519(94) 
0.3580(94) 
0.3051(40) 



0.2418(31) 
0.2448(31) 
0.1918(19) 
0.1975(19) 
0.1592(80) 
0.3494(82) 
0.3552(82) 
0.3059(42) 



TABLE XV: Heavy-quark mass dependence of radial and orbital energy splittings in lattice units, calculated with the «■* action 
on the L = 32 ensemble with ami ~ 0.004. 



3. Spin-dependent energy splittings 



ami = 0.005 



ami 



0.01 



ami 



0.02 



ami 



0.03 



r{lS)^rjb{lS) 

r{2S)-ribi2S) 

Xb2{lP) - XbiilP) 

XtijlP) - Xbo(lP) 

I'iP-hbilP) 

-2X60 (IP) - 3X61 (IP) + 5X62 (IP) 

-2x6o(lP) + 3x6i(lP)-Xw(lP) 

T(25')-r,f,(2S) 



T(IS') 



r?f,(lS') 
T(IS')- 



VbilS) 



-2xM,(l^) + 3x6i(l^)-X62(lP) 

y{2S)-r)b{2S) 

-2x6o(lP) + 3x6i(lP)-X62(lP) 



0.030216(78) 


0.03037(10) 


0.03087(23) 


0.03132(24 


0.0124(19) 


0.0135(22) 


0.0160(33) 


0.0150(27) 


0.01116(97) 


0.01185(76) 


0.0135(19) 


0.0115(24) 


0.01573(79) 


0.01540(73) 


0.0179(15) 


0.0150(20) 


0.00116(56) 


0.00139(46) 


0.0023(13) 


0.0035(16) 


0.0873(52) 


0.0900(43) 


0.103(10) 


0.087(13) 


0.0204(19) 


0.0189(15) 


0.0223(31) 


0.0183(48) 


0.411(63) 


0.446(73) 


0.52(11) 


0.480(87) 


1.48(14) 


1.60(13) 


1.38(19) 


1.71(44) 


0.61(11) 


0.72(13) 


0.72(18) 


0.82(26) 



TABLE XVI: Spin-dependent energy splittings in lattice units, from the v"^ action on the L = 24 lattices, for arrib ~ 2.536. 
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arm = 0.005 


arm = 0.01 


ami = 0.02 


ami = 0.03 


T(lS')-r,f,(lS) 


0.025607(67) 


0.025781(84) 


0.02620(19) 


0.02657(19) 


r{2S)~vti2S) 


0.0090(16) 


0.0102(19) 


0.0120(27) 


0.0113(25) 


Xb2{lP) - XtiilP) 


0.00910(80) 


0.00974(69) 


0.0108(19) 


0.0093(22) 


XbiilP) - XboilP) 


0.01396(71) 


0.01380(64) 


0.0158(14) 


0.0134(17) 


l^P-ht{lP) 


0.00085(50) 


0.00099(46) 


0.0018(11) 


0.0026(15) 


-2X60 (IP) - 3X61 (IP) + 5X62 (IP) 


0.0734(45) 


0.0763(40) 


0.086(11) 


0.073(13) 


-2X60(1P) + 3x6i(lP) - Xb2{lP) 


0.0189(16) 


0.0178(14) 


0.0206(33) 


0.0174(40) 


r{2S)-Vbi2S) 

r{is)-vbiis) 


0.350(64) 


0.395(75) 


0.46(10) 


0.425(95) 


T(1S)- 776(15) 


1.36(12) 


1.45(12) 


1.27(20) 


1.52(35) 


-2X6()(1P) + 3x6i(lP) - X62(1P) 


T(25) - 776(25) 


0.476(96) 


0.57(12) 


0.58(16) 


0.65(21) 


-2x6o(lP) + 3X6l(lP) - X62(1P) 



TABLE XVII: Spin-dependent energy splittings in lattice units, from the v^ action on the L = 24 lattices, for ami, ~ 2.536. 



ami =0.004 



ami = 0.006 



ami = 0.008 



T(15)- 776(15) 
T(25)- 776(25) 

X62(lP)-X6l(lP) 

X6i(lP)-X6o(lP) 

l^P-hbilP) 

-2X60 (IP) - 3X61 (IP) + 5X62 (IP) 
-2X60(lP)+3X6l(lP)-X62(lP) 

T(25)- 776(25) 



T(15) 



r76(15) 

T(15) - 776(15) 



-2x6o(lP) + 3x6i(lP)-X62(lP) 

T(25)- 776(25) 

-2x6o(lP) + 3x6i(lP)-X62(lP) 



0.024441(74) 


0.024408(74) 


0.02455(11) 


0.0114(12) 


0.0109(19) 


0.0127(13) 


0.00856(93) 


0.00910(98) 


0.0087(11) 


0.01268(78) 


0.01286(85) 


0.01321(78) 


0.00030(43) 


0.00104(56) 


0.00049(63) 


0.0681(54) 


0.0714(55) 


0.0696(57) 


0.0168(15) 


0.0168(19) 


0.0177(18) 


0.465(51) 


0.448(79) 


0.520(53) 


1.45(13) 


1.46(17) 


1.38(14) 


0.675(95) 


0.65(14) 


0.72(10) 



TABLE XVIII: Spin-dependent energy sphttings in lattice units, from the «■* action on the i = 32 lattices, for a7n6 = 1.87. 



ami = 0.004 



ami = 0.006 



ami = 0.008 



T(15)- 776(15) 
T(25)- 776(25) 
X62(1P)-Xm(1P) 
X6i(lP)-X6o(lP) 
l^P-hbilP) 

-2x6o(lP) - 3X6l(lP) + 5X62(1P) 
-2X60(lP)+3X6l(lP)-X62(lP) 

T(25)- 776(25) 
T(15)- 776(15) 

T(15)- 776(15) 

-2x6o(lP) + 3x6i(lP)-Xw(lP) 

T(25) - 776(25) 

-2x6o(lP) + 3x6i(lP)-X62(lP) 



0.020215(69) 

0.0083(12) 

0.00677(85) 

0.01122(69) 

0.00018(42) 

0.0562(49) 

0.0157(14) 

0.410(61) 

1.29(11) 

0.528(92) 



0.020170(59) 

0.0082(14) 

0.00735(79) 

0.01116(70) 

0.00075(52) 

0.0591(44) 

0.0150(16) 

0.404(70) 

1.34(14) 

0.54(11) 



0.020292(96) 

0.0096(12) 

0.00710(95) 

0.01141(68) 

0.00024(50) 

0.0582(53) 

0.0157(16) 

0.476(61) 

1.29(13) 

0.613(99) 



TABLE XIX: Spin-dependent energy splittings in lattice units, from the v^ action on the L — 32 lattices, for a7n6 = 1.87. 
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4. Dependence of spin splittings on the couplings C3 and C4 





C3 = 0.8 


C3 = 1.2 


C4 = 0.8 


C4 = 1.2 


r{lS)-rjtilS) 


0.98016(18) 


1.02148(19) 


0.67151(53) 


1.3808(12) 


r{2S)~vti2S) 


0.983(87) 


1.025(91) 


0.68(10) 


1.35(14) 


Xt2{lP) - XbiilP) 


0.832(49) 


1.162(54) 


1.020(62) 


0.954(61) 


XbiilP) ~ Xto{lP) 


0.935(38) 


1.064(35) 


0.786(34) 


1.239(48) 


r^p-ht{ip) 


0.95(37) 


1.01(43) 


0.67(36) 


1.32(63) 


-2XM)(1P) - 3X61 (IP) + 5X62 (IP) 


0.871(29) 


1.129(31) 


0.936(32) 


1.059(39) 


-2X60 (IP) + 3X61 (IP) - X62(1P) 


0.991(84) 


1.008(76) 


0.658(67) 


1.40(11) 


T(2S')-r,6(2S) 
T(lS')-r,,(15) 


1.003(89) 


1.003(89) 


1.02(15) 


0.98(10) 


T(1S)-»?6(15) 


0.989(83) 


1.013(78) 


1.02(10) 


0.989(77) 


-2X60(1P)+3X61(1P)-X62(1P) 


T(2S) - 776(25') 


0.99(16) 


1.02(15) 


1.05(23) 


0.98(16) 


-2x6o(lP) + 3x6i(lP)-X62(lP) 



TABLE XX: Dependence of the spin splittings, calculated with the «■* action, on the couplings C3 and C4. Shown is the ratio 
of the splitting with either C3 7^ 1 or C4 7^ 1 to the splitting with all Ci = 1, calculated using bootstrap. The data are for the 
L = 24 ensemble with ami = 0.005 and anib = 2.536. 





C3 = 0.8 


C3 = 1.2 


C4 = 0.8 


C4 = 1.2 


T(lS')-r,6(lS) 


0.97788(17) 


1.02411(20) 


0.64656(47) 


1.4180(11) 


r{2S)~r,t{2S) 


0.98(13) 


1.03(13) 


0.63(12) 


1.44(19) 


X62(lP)-X6l(lP) 


0.795(53) 


1.205(58) 


1.016(55) 


0.961(68) 


X6l(lP)-X60(lP) 


0.924(34) 


1.071(30) 


0.765(34) 


1.265(53) 


l^P~ht{lP) 


0.92(42) 


1.02(49) 


0.66(41) 


1.36(66) 


-2X60 (IP) - 3X61 (IP) + 5X62 (IP) 


0.845(28) 


1.154(32) 


0.920(29) 


1.077(40) 


-2x6o(lP) + 3x6i(lP)-X62(lP) 


0.987(71) 


1.006(62) 


0.641(59) 


1.41(11) 


r{2S)-nb{2S) 

riis)-vb{is) 


1.00(13) 


1.00(13) 


0.97(19) 


1.01(14) 


T(15)- 776(15) 


0.991(75) 


1.018(62) 


1.008(95) 


1.002(74) 


-2x6o(lP) + 3x6i(lP)-X62(lP) 


T(2S)-77b(25) 


1.00(21) 


1.04(19) 


0.96(26) 


1.02(20) 


-2X60(lP)+3X6l(lP)-X62(lP) 



TABLE XXL Dependence of the spin splittings, calculated with v^ action, on the couplings C3 and C4. 
the splitting with either C3 7^ 1 or C4 7^ 1 to the splitting with all Ci = 1, calculated using bootstrap. 
L = 24 ensemble with 0771; = 0.005 and ami, = 2.536. 



Shown is the ratio of 
The data are for the 



5. Heavy-quark mass dependence of spin splittings 



24 



arUb 



2.3 



anib 



2.536 



anib = 2.7 



T(2S')-r;f,(2S) 
Xb2{lP) - XbiilP) 
XtiilP) - XtoilP) 
13p~hb{lP) 

-2xbo{lP) - 3xm(1P) + 5X62(1P) 
-2xbo(lP) + 3x6i(lP)-XM(lP) 



0.0138(21) 

0.0122(11) 

0.01734(86) 

0.00120(62) 

0.0956(62) 

0.0225(22) 



0.0124(19) 

0.01116(97) 

0.01573(79) 

0.00116(56) 

0.0873(52) 

0.0204(19) 



0.0116(18) 

0.01060(86) 

0.01469(71) 

0.00120(54) 

0.0824(46) 

0.0189(18) 



TABLE XXII: Heavy-quark mass dependence of spin splittings in lattice unit s, from tlie u* action on tlie L ■ 
ami = 0.005 (Tlie results for tiie T(IS') - 77b(15') splitting are given in Table |XXlil I . 



24 ensemble with 



arrib = 2.3 



arrib — 2.536 



arrib — 2.7 



Fit result 



ami = 0.005 
ami = 0.01 
ami = 0.02 
am,i = 0.03 



0.032685(90) 
0.03291(11) 
0.03344(25) 
0.03403(26) 



0.030216(78) 
0.03037(10) 
0.03087(23) 
0.03132(24) 



0.028727(75) 
0.02885(10) 
0.02933(24) 
0.02968(23) 



A = 0.0615(10), B = 0.00597(41) 
A = 0.06295(90), B = 0.00554(34) 
A = 0.0636(28), 5 = 0.0058(11) 
A = 0.0675(20), B = 0.00470(79) 



TABLE XXIII: Heavy-quark mass dependence of the IS* hyperfine splitting in lattice units, from the v'^ action on the L = 24 



ensembles. In the last column of the table, the results of correlated fits using the functional form ai5x(i5) 
A/{amb) + B are shown. 



aE, 



'7b (IS) 



anib — 2.3 



anib = 2.536 



am,b = 2.7 



T(25')-r7f,(2S) 

Xb2{lP) - XbiilP) 

XbijlP) - XbojlP) 

13p~hb{lP) 

-2xw(lP) - 3x6i(lP) + 5X62(1P) 

-2xbo{lP) + 3xbi{lP)-Xb2{lP) 



0.0097(18) 

0.00964(87) 

0.01517(80) 

0.00083(54) 

0.0786(47) 

0.0208(19) 



0.0090(16) 

0.00910(80) 

0.01396(71) 

0.00085(50) 

0.0734(45) 

0.0189(16) 



0.0085(15) 

0.00877(83) 

0.01320(64) 

0.00093(52) 

0.0703(47) 

0.0177(15) 



TABLE XXIV: Heavy-quark mass dependence of spin splittings in lattice units, from the v^ action on the L = 24 ensemble 
with ami = 0.005 (The results for the T(IS') - rjbilS) splitting are given in Table|XXV"l. 



arrib = 2.3 



arrib = 2.536 



arrib = 2.7 



Fit result 



ami = 0.005 
ami = 0.01 
ami = 0.02 
arr%i = 0.03 



0.027319(70) 
0.027548(91) 
0.02800(19) 
0.02851(21) 



0.025607(67) 
0.025781(84) 
0.02620(19) 
0.02657(19) 



0.024567(66) 
0.024712(85) 
0.02513(19) 
0.02540(19) 



A = 0.0427(10), B = 0.00876(40) 
A = 0.04398(75), B = 0.00843(29) 
A = 0.0447(25), 5 = 0.0086(10) 
A = 0.0482(16), B = 0.00753(65) 



TABLE XXV: Heavy-quark mass dependence of the 15* hyperfine splitting in lattice units, from the v action on the L = 24 
ensembles. In the last column of the table, the results of correlated fits using the functional form aE-r(is) ^ 0^776(15) = 
A/{amb) + B are shown. 



25 



T(25')-r?f,(2S) 
Xb2{lP) - XbiilP) 
XbijlP) - XbojlP) 
l^P-htilP) 

-2xbo{lP) - 3xbi(lP) + 5xb2{lP) 
-2xbo{lP) + SXbiilP) - Xb2{lP) 



arnt = 1.75 

0.0122(13) 

0.0091(10) 

0.01352(84) 

0.00028(49) 

0.0725(60) 

0.0179(18) 



arub = 1.87 

0.0114(12) 

0.00856(93) 

0.01268(78) 

0.00030(43) 

0.0681(54) 

0.0168(15) 



arub ~ 2.05 

0.0108(12) 

0.00804(76) 

0.01172(70) 

0.00033(39) 

0.0636(44) 

0.0154(13) 



TABLE XXVI: Heavy-quark mass dependence of spin splittings in lattice units, from the v^ action on the L 
with ami = 0.004 (The results for the T(15') - 776(15) splitting are given in Table |XXV"n I . 



32 ensemble 



arub = 1.75 arrib — 1.87 arUb — 2.05 Fit result 

ami = 0.004 0.025679(82) 0.024441(74) 0.022850(76) A = 0.03380(58), B = 0.00636(28) 

ami = 0.006 0.025669(78) 0.024408(74) 0.022787(70) A = 0.03443(54), B = 0.00598(27) 

am,i = 0.008 0.02582(12) 0.02455(11) 0.02291(11) A = 0.0348(11), B = 0.00592(57) 

TABLE XXVII; Heavy-quark mass dependence of the IS hyperfine splitting in lattice units, from the n* action on the L = 32 

ensembles. In the last column of the table, the results of correlated fits using the functional form aET(is) — 'i^)7b(is) = 
A/{amb) + B are shown. 



arUb 



1.75 



arUb 



1.87 



arub 



2.05 



T(25')-r;f,(2S) 
Xb2{lP) ~ Xbi{lP) 
XbijlP) - XbojlP) 
13P~hb{lP) 

~2xbo{lP) - SxbiilP) + 5xb2{lP) 
-2xbQ{lP) + 3xbi{lP)-Xb2{lP) 



0.0088(13) 

0.00700(89) 

0.01187(73) 

0.00020(42) 

0.0587(51) 

0.0168(15) 



0.0083(12) 

0.00677(85) 

0.01122(69) 

0.00018(42) 

0.0562(49) 

0.0157(14) 



0.0077(13) 

0.00643(72) 

0.01040(61) 

0.00020(32) 

0.0530(41) 

0.0144(12) 



TABLE XXVIII: Heavy-quark mass dependence of spin splittings in lattice units, from t he v action on the L 
with ami = 0.004 (The results for the T(15) - Vbi'i^S) splitting are given in Table|XXIXl. 



32 ensemble 



amt 



1.75 



amb 



1.87 



arUb 



2.05 



Fit result 



ami = 0.004 0.021117(74) 0.020215(69) 0.019088(62) A = 0.02424(44), B = 0.00726(22) 

ami = 0.006 0.021083(63) 0.020170(59) 0.019029(56) A = 0.02450(56), B = 0.00707(27) 

ami = 0.008 0.021222(92) 0.020292(96) 0.019126(91) A = 0.0251(11), B = 0.00685(58) 



TABLE XXIX: Heavy-quark mass dependence of the IS hyperfine splitting in lattice units, from the v action on the L = 32 



ensembles. In the last column of the table, the results of correlated fits using the functional form a_Ex(is) 
A/{arab) + B are shown. 



a£, 



VbiiS) 
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FIG. 9: Heavy-quark mass dependence of the spin splittings on the L — 24, ami = 0.005 and L — 32, ami = 0.004 ensembles. 
Results are shown for both the u* and the v^ actions. The dashed lines and gray error bands are correlated fits using the 
function A/{amb). The data for the 15 hyperfine splittings are incompatible with this form, and additional fits using the 
function A/{amt) + B axe shown, which describe the data very well. 
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Appendix D: Gluon discretization errors 

The Iwasaki giuon action used in this work belongs to a class of actions with the form 



Sg[U] 



/3 



E 



(1 - 8ci) V PpU^,, + ci V R[U]., 



tJ.<v 



fl^U 



(Dl) 



where P[U]x,^iy and R[U]x.fj.i, are the real part of the trace of the 1x1 plaquette and 1x2 rectangle terms, respectively 
(the coefficient ci in ( |D1[ ) should not be confused with the one in the NRQCD action Q). The Iwasaki action uses 
ci = —0.331, derived from a renormalization-group transformation [28 -31?] . Note that tree-level order-a^ improvement 
would require ci = —1/12, corresponding to the tree-level Liischer-Weisz action (IHlUnj)- However, nonperturbatively 
and at coarse lattice spacings, the Iwasaki action has been shown to yield reduced lattice artifacts compared to the 
tree-level Liischer-Weisz action [SDJ [5T] . 

In this appendix, gluonic discretization errors in bottomonium energy splittings will be investigated. In Sec. |D 1[ 
the shifts in radial and orbital energy splittings are studied using tree-level perturbation theory for four different 
choices of ci. The tree-level energy shift for the simple plaquette action (ci = 0) has previously been estimated at 
order a^ in [S21[S3]. In the following, a new analysis based on a lattice potential model is presented. This analysis 
does not make use of an expansion in powers of a, which would not be appropriate for the Iwasaki action. 



Then, to go beyond tree-level, in Sec. D 2 nonperturbative bottomonium results obtained from the RBC/UKQCD 
ensembles (using the Iwasaki gluon action) and from the MILC ensembles '13' (using the tadpole-improved one-loop 
Liischer-Weisz action |54| ) are compared. This comparison also includes all the bottomonium spin splittings considered 
in this paper, and leads to nonperturbative estimates of gluonic discretization errors for them. 



1. Lattice potential model using tree- level perturbation theory 

a. The model 

The discretization errors caused by the gluon action in radial and orbital bottomonium energy splittings can be 
estimated using a potential model on a three-dimensional cubic lattice with Hamiltonian 



H 



nib 



V, 



(D2) 



where A is a lattice Laplace operator and V{r) is the static quark-antiquark potential derived from the lattice gluon 
action in use. In the model employed here, V{r) is taken to be of the form 



V{r) = V''''-'"{r) + K\rl 



(D3) 



where V ' {r) is the tree-level lattice potential that is obtained from the tree-level lattice gluon propagator G!^l{q) 
as follows: 



yIat,U(^) 



d3 



i^ 4l<-(2^)'' 



^ iq-r ^lat 



G;,1f(q,0). 



(D4) 



For a — >■ 0, the potential y''^''°(r) approaches the continuum Coulomb potential V^{r) = — (4/3)as/|T| (with ag — 
g^ /{An)). The linear term K\r\ in (D3) is added to describe the nonperturbative long-distance behavior of the quark- 
antiquark potential. In Ref. [39] a similar model on a cubic lattice was considered, with a discrete Laplacian but 
with the continuum form of the potential. In the following, the same parameters as in [5S] are used: ^/K = 468 MeV, 
as = 0.24, and rrib — 4.676 GeV. Note that in ^39^ the coordinate system was chosen such that the origin r = 
was at the center of an elementary cube, in order to avoid the singularity of the continuum Coulomb potential. In 



contrast, here the point r = is a regular lattice point, and the lattice potential (D4| is finite at that point (one has 
]/iat,o^Q-j ^ 1/a). In fact, the dominant g luon discretization errors arise at and near the point r = 0. 



The propagator Gjf*(g) for the action (Dl) can be found in 55]. The 0-0-component at go = has the simple form 



4o(g,o) = 



{lYT. 



{-. 



Clfl- 



(D^E 



{-. 



(D5) 
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The Iwasaki gluon action p5"W has ci = —0.331. Results will also be given for Ci = (the simple plaquette action), 
ci = -1/12 (the tree-level Liischer-Weisz action j48l|49]), and ci = -1.40686 (the DBW2 action [SI IS])- For these 
choices of Ci and for all points r with \ri/a\ < 60 the integral (D4) was computed numerically. 
For the Laplace operator A in ( D2 ) , three different discretizations are considered: 



fE 



=iv;v7' 



^-^J = l J 3 

=1 J J 



T 



(«Vi2)E-=i[v+v-]', 

(«V12) E -=1 [V+ V-] ' + (aV90) EU M^J] ' ' 



unimproved, 
O (a^ )-improved, 
C'(a^)-improved, 



(D6) 



with V^ip{r) = [ijj{r + aej) — ijj{r)] /a and V ■ ijj{r) ~ [ipir) — ijj{r — aSj)] /a. As shown by the results in the next 
section, when the C'(a'*)-improved Laplacian is used, the discretization errors associated with A are in most cases 
much smaller than the gluonic discretization errors associated with ( D4 1 . 

The low-lying eigenvalues and eigenfunctions of the Hamiltonian ( D2 ) were computed numerically for lattices with 



a physical side length of 2.7 fm and lattice spacings in the range from 0.0223 fm to 0.208 fm. As in IW' , only one 
octant of the lattice was simulated. S'-wave states (Ai representation) and Z?-wave states {E representation) were 
obtained by using periodic boundary conditions in all three lattice directions; P-wave states (Ti representation) were 
obtained by using antiperiodic boundary conditions in the ra-direction and periodic boundary conditions in the ri- 
and r2-directions. 



Results 



Figure 10 shows the deviations in the IS" and IP energies from their continuum values as a function of a^, for 
the three different levels of Symanzik improvement in the lattice Laplacian defined in Eq. ( D6 ) . As can be seen in 



the figure, the difference between the results from the 0{a^)- and C'(a^)-improved Laplace operators is small, much 
smaller than the difference between the results from the unimproved and C'(a^)-improved Laplace operators. This 
indicates that the remaining discretization errors in the C'(a^)-improved Laplacian are negligible. However, for the 
S'-wave state a significant shift in the energy from its continuum value remains at finite lattice spacing. This error can 
be interpr eted as the tree-level gluonic contribution to the discretization errors, stemming from the use of the lattice 
potential (D4|. As expected, at small lattice spacings, the C'(a^)-improved Liischer-Weisz action shows significantly 



smaller tree-level discretization errors than the other actions. 

For the P-wave states, the gluonic discretization errors are much smaller than for the S wave states. This is expected 
because the dominant correction in the potential arises at the origin, where the wave function vanishes for all states 
other than the S-wave states. Only the DBW2 action leads to significant tree-level gluon discretization errors in the 
IP energy. 



Examples of IS" and IP eigenfunctions of the Hamiltonian (D2) with the C'(a^)-improved Laplacian are shown in 
Fig. 11 The large negative coefficient ci — —1.40686 of the DBW2 action leads to a visible distortion of the 15 wave 
function compared to the tree-level Liischer-Weisz action with ci ~ —1/12, while the IP state is only weakly affected 
by the choice of ci. The broadening of the IS* wave function is expected because a negative coefficient ci shifts the 
potential at short distances upwards. In particular, at r = 0, one has 



3a 



^lat,U(o) 



-0.25273 , 
-0.21903, 
-0.16437, 
-0.09384 , 



ci=0, 
ci = -l/12, 
ci = -0.331, 
ci = -1.40686. 



(D7) 



However, note that the broadening of the lattice 15 wave function caused by a negative value of ci does not necessarily 



mean that the hyperfine splitting is reduced. This will be discussed further in Sec. D 2 



The results for the shifts in the energies of the 15, 25, 35, IP, 2P and ID states, obtained with the C'(a'*)-improved 
Laplacian and the four different gluon actions, are summarized in Fig. [121 This figure also shows the 25 — 15 and 
IP — 15 splittings, demonstrating that the 25 — 15 splitting has smaller tree-level gluon discretization errors and is 
therefore better suited for setting the lattice scale. For the Iwasaki action, the gluonic tree-level discretization errors 
in the 25 — 15 splitting are found to be about 2.6% at a = 0.11 fm and 1.6% at a = 0.08 fm, respectively. Note that 
the 2P — IP splitting is nearly free of gluonic discretization errors and therefore appears to be a good alternative 
choice for the scale setting. However, in the actual lattice QCD calculation the 2P — IP splitting has much larger 
statistical errors than the 25 — 15 splitting (see Sec. C 1 1. 
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FIG. 10: Shift in the lattice potential model IS and IP energies as a function of the lattice spacing, for different levels of 
Symanzik improvement in the Laplace operator. 
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FIG. 11: Lattice potential model wave functions '^{ri,r2,r3) at ri = for the IS and IP states. Data are shown for the 
tree-level Liischer-Weisz gluon action (ci = —1/12) and the DBW2 gluon action (ci — —1.40686), at two different lattice 
spacings, using the C'(a'')-improved Laplacian in all cases. 
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FIG. 12; Shift in the lattice potential model energy levels and sphttings as a function of the lattice spacing. Ah data shown in 
this figure were generated with the 0(a'')-improved Laplacian, so that the shifts are dominated by gluon discretization errors. 
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FIG. 13: Shift in the lattice potential model energy splittings, rescaled using the 25 — 15 splitting, as a function of the lattice 
spacing. AH data shown in this figure were generated with the 0(a'')-iinproved Laplacian, so that the shifts are dominated by 
gluon discretization errors. 

Figure [T3| shows the gluonic discretization errors in the 3S — IS, IP — IS, 2P — IP, and ID — IS sphttings for the 
case that they are calculated using the 25" — IS" splitting to set the scale. In the 2P — IP splitting, previously nearly 
free of gluonic discretization errors, this process introduces new gluon errors. However, in the other splittings shown 
in Fig. [13] the scale setting with the 25' — IS* splitting leads to a partial cancellation of gluonic discretization errors. 
The results for the relative errors obtained with the Iwasaki action at a = 0.11 fm and a = 0.08 fm are given in Table 

Ixxxl 



error (a = 0.11 fm) 



error (a = 0.08 fm) 



35-15 


0.6% 


IP- IS 


-3.4% 


2P- IP 


2.6% 


ID -15 


-1.2% 



0.4% 
-2.8% 

1.6% 
-1.1% 



TABLE XXX: Estimates of tree-level gluon discretization errors in radial and orbital bottomonium energy splittings computed 
with the Iwasaki action and using the 25 — 15 splitting to set the scale. A negative sign indicates a negative deviation from 
the continuum value. 



2. Results from MILC ensembles and gluon discretization errors in spin splittings 



In order to study the influence of the gauge action on the bottomonium energy splittings nonperturbatively, the 
calculations presented in the main part of this paper were repeated on two ensembles of lattice gauge fields generated by 
the MILC collaboration [13]. These ensembles make use of the tadpole-improved one- loop Liischer-Weisz action [54] for 
the gluons, which is based on order-a^ Symanzik improvement rather than renormalization-group improvement. The 
action includes the plaquette and rectangle terms, and in addition a third term ("parallelogram"). Their coefficients 
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L^xT 


/3 


ami 


arus 


arub 


UOL 


^conf 


a (fm) 


mf ^S' (GeV) 


(phys.) 


24^ X 64 
28^ X 96 


6.76 
7.09 


0.005 
0.0062 


0.05 
0.031 


2.3, 2.64, 2.7 
1.75, 1.86, 2.05 


0.8362 
0.8541 


2099 
1910 


0.1198(11) 
0.08613(83) 


0.460 
0.416 


2.664(24) 
1.858(19) 



TABLE XXXI: Parameters for the calculations on the MILC ensembles. The bare gauge couplings are given as P — IQ/g^ . For 
the MILC ensembles, there are taste splittings between the different pions [T3], and the root-mean-square masses taken from 
[58n59j are given. 



/3pij Pit, and /3pg were computed using one- loop perturbation theory, but without the effects of sea quarks. These 
effects were later calculated and found to be significant [57] ■ Therefore, on the (2 + l)-flavor MILC ensembles, the 
gluon action is expected to have 0{asa^) errors. 

The parameters of the MILC ensembles used here are given in Table |XXXI[ The values for the lattice spacing 
and arrv^ ^^ were computed using exactly the same methods as for the RBC/UKQCD ensembles, to minimize any 
possible bias. The sea quarks in the MILC ensembles are implemented with the rooted staggered AsqTad action 
[50HB2]- This leads to an effective averaging over multiple tastes of sea pions [T3], and therefore the appropriate 
pion mass to consider for bottomonium is the root-mean square (RMS) pion mass. To facilitate the comparison, the 
results from the RBC/UKQCD ensembles were therefore interpolated/extrapolated to match the RMS pion masses 
of the MILC ensembles. The lattice spacings of the coarse and fine MILC ensembles also match the lattice spacings 
of the corresponding RBC/UKQCD ensembles. As discussed in Sec. |IIIB2 the lattice spacing of the RBC/UKQCD 



ensembles changes slightly when the sea quark mass is changed, because the bare gauge coupling is kept constant. 
It turns out that this shift makes the agreement of the lattice spacings even better after interpolation to match the 
MILC pion masses. At the matching points, any significant difference between the results from the RBC/UKQCD 
and MILC ensembles would indicate different systematics associated with the gluon and sea quark actions. In the 
following it is assumed that the effect of changing the sea quark action is negligible for bottomonium. 



The radial and orbital energy splittings are compared in Table XXXII As can be seen there, with the exception 
of the statistically most precise IP — IS splitting, all results from the MILC ensembles agree with those from the 
RBC/UKQCD ensembles within the statistical errors. At the coarse lattice spacing, the IP — IS splitting from the 
MILC ensemble is found to be about 9 MeV (1.4 standard deviations) higher than that from the RBC/UKQCD 
ensemble. At the fine lattice spacing, the IP — IS splittings from the MILC and RBC/UKQCD ensembles fully agree 
with each other within the statistical error of 7 MeV. In contrast, the tree-level estimates in Fig. [T3] would suggest a 
difference between the splittings from the Iwasaki and tree-level Liischer-Weisz actions of 27 MeV at a = 0.12 fm and 
18 MeV at a = 0.09 fm. Note however that at these lattice spacings, the tadpole-improved one-loop Liischer-Weisz 
action has a value of (3rt/l3pi that is not as far away from the Iwasaki action as in the tree- level case [54]. Thus, 
nonperturbatively the errors caused by the Iwasaki action are likely to be smaller than the tree-level estimates. 



a Ri 0.12 fm, m^ ^ 460 MeV 



a ^ 0.09 fm, m^ ^ 416 MeV 



RBC 



MILC 



RBC 



MILC 



T(35') - T(IS') 
]AP-T(1S') 

Wp-Ts 

23p_l3p 

23P-T(1S') 

Wp-Ts 

T2{lD)-r{lS) 



0.916(28) 

0.4426(41) 

0.4556(41) 

0.3932(92) 

0.836(11) 

0.849(11) 

0.7130(70) 



0.915(30) 

0.4510(45) 

0.4645(46) 

0.385(12) 

0.831(12) 

0.845(13) 

0.7143(70) 



0.926(23) 

0.4349(45) 

0.4490(46) 

0.396(13) 

0.829(14) 

0.843(14) 

0.7091(81) 



0.934(31) 

0.4367(58) 

0.4515(60) 

0.386(11) 

0.823(15) 

0.837(15) 

0.713(12) 



TABLE XXXII: Comparison of results from the RBC/UKQCD and MILC ensembles: radial and orbital energy splittings in 
GeV, computed with the w^ action. 



XXXIII 



The spin splittings from t he RBC/ UKQCD and MILC ensembles are compared in Table 
NRQCD action) and Table |XXXLV| (for the order-w^ NRQCD action). As can be seen there 



(for the order-w^ 



with the exception 

of the directly calculated 15 hyperfine splitting, all results from the MILC ensembles are in agreement with the 
corresponding results from the RBC/UKQCD ensembles within the statistical errors. At the fine lattice spacing, the 
15' hyperfine splitting from the MILC ensemble is found to be about 5% (1.9 standard deviations) higher than the 
15 hyperfine splitting from the RBC/UKQCD ensemble. 
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a ^ 0.12 fm, 77i,r 


^ 460 MeV 


a ^ 0.09 fm, m„ 


^ 416 MeV 




RBC 


MILC 


RBC 


MILC 


T(15)- 776(15) 


52.29(99) 


53.79(95) 


56.2(1.1) 


59.3(1.2) 


r{2S)--nbi2S) 


23.3(2.1) 


25.0(2.2) 


28.1(2.1) 


31.6(3.2) 


Xb2{lP) - XbiilP) 


20.3(1.1) 


21.5(1.3) 


20.4(1.5) 


21.8(1.8) 


XbiilP) ~ Xbo{lP) 


27.1(1.0) 


27.3(1.1) 


29.8(1.3) 


29.8(1.6) 


l^p-htilp) 


2.94(66) 


2.90(70) 


1.93(84) 


2.1(1.2) 


-2xw(lP) - 3x61 (IP) + 5x62 (IP) 


155.3(6.4) 


162.1(7.5) 


161.6(8.8) 


168.7(9.6) 


-2X60(lP)+3X6l(lP)-X62(lP) 


33.9(2.1) 


33.1(2.3) 


39.4(2.7) 


38.0(3.8) 


T{2S)-Vb{2S) 


0.447(39) 


0.466(41) 


0.501(38) 


0.532(52) 


r{is)~vb{is) 


1.517(88) 


1.63(11) 


1.416(98) 


1.56(15) 


-2x6o(lP) + 3x6i(lP)-X62(lP) 


T(2S)- 7,6(25) 


0.692(75) 


0.755(83) 


0.720(73) 


0.83(12) 


-2x6o(lP)+3X6i(lP)-X62(lP) 



TABLE XXXIII: Comparison of results from the RBC/UKQCD and MILC ensembles: spin splittings, computed with the v'^ 
action. All results in MeV, except for the dimensionlcss ratios. 





a « 0.12 fm, TTivr 


« 460 MeV 


a « 0.09 fm, 7n^ 


^ 416 MeV 




RBC 


MILC 


RBC 


MILC 


T(15)- 776(15) 


44.39(84) 


45.64(81) 


46.49(91) 


48.64(96) 


T(25)- 776(25) 


17.3(1.8) 


18.8(1.9) 


20.9(1.9) 


22.3(2.4) 


X62(lP)-X6l(lP) 


16.55(99) 


17.7(1.2) 


16.5(1.3) 


17.5(1.7) 


X6l(lP)-X60(lP) 


24.06(91) 


24.3(1.0) 


25.9(1.1) 


26.5(1.5) 


l^P^htilP) 


2.18(60) 


2.13(62) 


1.27(74) 


1.21(96) 


-2X60 (IP) - 3X61 (IP) + 5X62 (IP) 


130.6(5.9) 


137.0(6.9) 


134.3(7.8) 


140.9(8.8) 


-2X60(lP)+3X6l(lP)-X62(lP) 


31.5(1.9) 


30.9(2.2) 


35.6(2.4) 


35.7(3.7) 


T(25)- 776(25) 
T(15)- 776(15) 


0.392(40) 


0.413(42) 


0.450(41) 


0.458(47) 


T(15) - 7,6(15) 


1.396(81) 


1.48(10) 


1.304(88) 


1.36(14) 


-2X60(1P)+3X61(1P)-X62(1P) 


T(25)- 776(25) 


0.552(66) 


0.607(75) 


0.593(67) 


0.626(90) 


-2X60(1P)+3X61(1P)-X62(1P) 



TABLE XXXIV: Comparison of results from the RBC/UKQCD and MILC ensembles: spin splittings, computed with the 
action. All results in MeV, except for the dimensionless ratios. 



Recall from Sec. |D 1 b| that a negative coefficient ci in the gluon action leads to a broadening of the lattice IS" 
wave function. In the continuum, the leading-order hyperfine splitting is proportional to |'0(O)p, and one might 
therefore expect naively that a negative coefficient ci reduces the hyperfine splitting [63 . If this picture was correct, 
for example the simple plaquette action (ci — 0) would give a significantly larger hyperfine splitting than the Iwasaki 
action (ci = —0.331). In Ref. [22], the authors compared their results for the IS hyperfine splitting, computed using 
the Iwasaki action, to the results from [3T] that used the same NRQCD action and the same number of sea quark 
flavors, but the plaquette gluo n action. The lattice spacing was a ~ 0.10 fm in both cases. At a = 0.10 fm, the lattice 



potential model from Sec. 



Dl 



gives 



|^(ci=O)(0)|2/|^, 



(ci=-0.331) 



(O)l^ 



L5. 



In contrast, the results for the IS" hyperfine 
splitting from the two groups were in agreement within the statistical error of about 5%. Clearly, the simple continuum 
picture for the hyperfine splitting does not apply here. On the lattice, the spin-dependent potential responsible for 
the S'-wave hyperfine splitting will have non-zero values also at r 7^ 0. It is expected to be a complicated function 
that depends both on the lattice gluon propagator and on the discretization of the chromomagnetic field strength in 
the NRQCD action. 

Given the comparison of results for the 15* hyperfine splitting from the Iwasaki action with results from two other 
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gluon actions as discussed above, it seems reasonable to assume a gluonic discretization error of 5% for the hyperfine 
splittings calculated with the Iwasaki action at a w 0.08 fm. The same error estimate is used for the ratios of hyperfine 
and tensor splittings. For the ratio of the 25* and 15' hyperfine splittings, a partial cancellation of gluonic discretization 
errors is expected (as in the 2S — IS splitting in Fig. 12), and therefore a 2.5% gluon error is estimated for this ratio 
at a « 0.08 fm. 

Finally, recall from Fig. [o] that the IS" hyperfine splitting computed on the RBC/UKQCD ensembles shows an 
arrifj-dependence that is slightly different from the simple proportionality to l/{ami,) seen in the other spin splittings. 
In Fig.[T4] the 15 hyperfine splitting in lattice units computed on the coarse and fine MILC ensembles is plotted as a 
fimction of l/{ami,). As can be seen there, the behavior is very similar to that found on the RBC/UKQCD ensembles. 
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FIG. 14: Heavy-quark mass dependence of the 15* hyperfine splittings on the MILC ensembles (see Fig. lo] for the data from 
the RBC/UKQCD ensembles). 



[11 
[12: 

[13: 

[14 
[15 
[16 
[17 
[18 
[19 
[20; 

[21 
[22 
[23 
[24 



[1] B. Aubert et al. (BABAR Collaboration), Phys. Rev. Lett. 101, 071801 (2008) [Erratum-ibid. 102, 029901 (2009)] 

arXiv: 0807. 1086 [hep-ex]]. 
[2] B. Aubert et al. (BABAR Collaboration), Phys. Rev. Lett. 103, 161801 (2009) |arXiv :0903. 1124| [hep-ex]]. 
[3] G. Bonvicini et al. (CLEG Collaboration), Phys. Rev. D 81, 031104 (2010) arXiv : 090 9 . 5474| [hep-ex]]. 
[4] S. Recksiegel and Y. Sumino, Phys. Lett. B 578, 369 (2004) arXiv:hep-ph/0305178 . 
[5] B. A. Kniehl, A. A. Penin, A. Pineda, V. A. Smirnov and M. Steinhauser, Phys. Rev. Lett. 92, 242001 (2004) [Erratum-ibid. 

104, 199901 (2010)] |arXiv:hep-ph/0312086 . 
[6] A. A. Penin, arXiv: 0905. 4296 [hep-ph] . 
[7] F. Domingo, U. Ellwanger and M. A. Sanchis-Lozano, Phys. Rev. Lett. 103, 111802 (2009) arXiv: 0907. 0348 [hep-ph]]. 

A. X. EbKhadra, A. S. Kronfeld and P. B. Mackenzie, Phys. Rev. D 55, 3933 (1997) |a rXiv:hep-lat/960400 4 . 
[9] B. A. Thacker and G. P. Lepage, Phys. Rev. D 43, 196 (1991). ' 

[10] G. P . Lepage, L. Magnea, G. Nakhleh, U. Magnea and K. Hornbostel, Phys. Rev. D 46, 4052 (1992) 



arXiv : hep-lat/9205007 



X. Liao and T. Manke, Phys. Rev. D 65, 074508 (2002) arXiv:hep-lat/0111049 . 

T. W. Chiu, T. H . Hsieh, C. H. Huang and K. Ogawa (TWQGD Collaboration), Phys. Lett. B 651, 171 (2007) 

[arXiv: 0705 .27971 [hep-lat]]. 

A. Bazavov et al, Rev. Mod. Phys. 82, 1349 (2010) arXiv: 0903. 3598' [hep-lat]]. 

T. Burch et al., Phys. Rev. D 81, 034508 (2010) arXiv: 0912. 2701 [hep-lat ]_[. 



M. Li (RBC and UKQCD Collaborations), PoS LATTICE2008, 120 (2008) arXiv: 0810. 0040 [hep-lat]]. 

C. AUton et al. (RBC and UKQCD Collaborations), Phys. Rev. D 78, 114509 (2008) arXiv: 08 04. 0473 [hep-lat]]. 

E. Eichten, K. Gottfried, T. Kinoshita, K. D. Lane and T. M. Yan, Phys. Rev. D 21, 203 (1980). 

A. Gray et al., Phys. Rev. D 72, 094507 (2005) ^arXiv:hep-lat/0507013 . 

H. D. Trottier, Phys. Rev. D 55, 6844 (1997) (arXiv: hep-lat/9611026 . 

T. Manke, I. T . Drummond, R. R. Horgan and H. P. Shanahan (UKQCD Collaboration), Phys. Lett. B 408, 308 (1997) 



arXiv : hep-lat/9706003 



N. Eicker et al, Phys. Rev. D 57, 4080 (1998) arXiv : hep-lat 79769002 . 

T. Manke et al. (CP-PACS Collaboration), Phys. Rev. D 62, 114508 (2000) |arXiv:hep-lat/0005022 

S. Meinel, Phys. Rev. D 79, 094501 (2009) |arXiv: 0903. 3224 [hep-lat]]. 



C. AUton et al. (RBC and UKQCD Collaborations), Phys. Rev. D 76, 014504 (2007) [arXiv : hep-lat/07010l3 



36 



[25 
[26 
[27 
[28 
[29 
[30 
[31 
[32 
[33 
[34 
[35 

[36 

[37 
[38 
[39 
[40 

[41 
[42 
[43 
[44 
[45^ 

[46 

[47; 

[48 
[49 
[50 
[51 
[52 

[53 

[54; 

[55 
[56 
[57 



D. B. Kaplan, Phys. Lett. B 288, 342 (1992) | arXiv : hep-lat/9206013] . 
Y. Shamir, Nucl. Phys. B 406, 90 (1993) arXiv : hep-lat / 9303005, . 



V. Furman and Y. Shamir, Nucl. Phys. B 439, 54 (1995) | larXiv:hep-lat/9405004 

Y. Iwasaki, Report No. UTHEP-118 (1983) 

Y. Iwasaki and T. Yoshie, Phys. Lett. B 143, 449 (1984) 



Y. Iwasaki, K. Kanaya, T. Kaneko and T. Yoshie, Phys. Rev. D 56, 151 (1997) [arXiv:hep-lat/9610023 
Y. Aoki et ai, Phys. Rev. D 69, 074504 (2004) arXiv:hep-lat/0211023 . 

S. N. Syritsyn et al, PoS LATTICE2008, 169 (2008) arXiv: 0903. 3063 [hep-lat]]. 

G. P. Lepage and P. B. Mackenzie, Phys. Rev. D 48, 2250 (1993) |arXiv:hep-lat/9209022; 



N. H. Shakespeare and H. D. Trottier, Phys. Rev. D 58, 034502 (1998) arXi v : hep-lat /980 2038| . 

G. P. Lepage, B. Clark, C. T. H. Davies, K. Hornbostel, P. B. Mackenzie, C. Morningstar and H. Trottier, Nucl. Phys. 

Proc. Suppl. 106, 12 (2002) |arXiv:hep-lat/0110175 . 

C. Amsler et al. (Particle Data Group), Phys. Lett. B 667, 1 (2008). 



G. Bonvicini et al. (CLEO Collaboration), Phys. Rev. D 70, 032001 (2 004) ,arXiv : hep-ex/040402 1 



P. del Amo Sanchez et al. (BABAR Collaboration), a rXiv: 1004.0175 [hep-ex]. 

G. S. Bah and P. Boyle, Phys. Rev. D 59, 114504 (1999) 'arXiv: hep-lat /98091801. 

G. T. Bodwin, E. Braaten and G. P. Lepage, Phys. Rev. D 51, 1125 (1995) [Erratum-ibid. D 55, 5853 (1997)] 

arXiv : hep-ph/9407339 . 

E. FoUana et al. (HPQCD and UKQCD Collaborations), Phys. Rev. D 75, 054502 (200 7) |arXiv:hep-lat/0610092 . 

C. S. Kim, T. Lee and G. L. Wang, Phys. Lett. B 606, 323 (2005) arXi v : hep-ph/04 1 imt^. 

E. B. Gregory et al., Phys. Rev. Lett. 104, 022001 (2010) arXiv: 0909. 4462 [hep-lat ]]""^ 



D. Ebert, R. N. Faustov and V. O. Galkin, Phys. Rev. D 67, 014027 (2003) arXiv : hep -ph/0210381 . 

A. M. Badalian, B. L. G. Bakker and I. V. Danilkin, Phys. Rev. D 81, 071502 (2010) [Erratum-ibid. D 81, 099902 (2010)] 

larXiv: 0911. 4634 [hep-ph]]. 

S. F. Radford and W. W. Repko, arXiv: 0912. 2259 [hep-ph]. 

R. G. Edwards and B. Joo (SciDAC, LHPC and UKQCD CoUaborations), Nucl. Phys. Proc. Suppl. 140, 832 (2005) 



arXiv : hep-lat/0409003 



M. Liischer and P. Weisz, Commun. Math. Phys. 97, 59 (1985) [Erratum-ibid. 98, 433 (1985)]. 
M. Liischer and P. Weisz, Phys. Lett. B 158, 250 (1985). 



P. de Forcrand et al. (QCD-TARO Collaboration), Nucl. Phys. B 577, 263 (2000) [arXiY:hep-lat/9911033 

S. Necco, Nucl. Phys. B 683, 137 (2004) arXiv:hep-lat/0309017 . 

C. T. H. Davies, K. Hornbostel, G. P. Lepage, A. Lidsey, J. Shigemitsu and J. H. Sloan, Phys. Lett. B 345, 42 (1995) 

,arXiv : hep-ph/9408328 . 

C. T. H. Davies, K. Hornbostel, G. P. Le page, A. Lidsey, P. McCal lum, J. Shigemitsu and J. H. Sloan (UKQCD Collabo- 
ration), Phys. Rev. D 58, 054505 (1998) 'arXiv : hep-lat /9802024' . 

M. G. Alford, W. Dimm, G. P. Lepage, G. Hockney and P. B. Mackenzie, Phys. Lett. B 361, 87 (1995) 
arXiv : hep-lat/9507010 

P. Weisz and R. Wohlert, Nucl. Phys. B 236, 397 (1984) [Erratum-ibid. B 247, 544 (1984)]. 

T. Takaishi, Phys. Rev. D 54, 1050 (1996). 

Z. Hao, G. M. von Hippel, R. R. Horgan, Q. J. Mason and H. D. Trottier, Phys. Rev. D 76, 034507 (2007) [arXiv : 0705 ■ 4660 

[hep-lat]]. ^__^__ 

[58] A. Bazavov et al. (MILC CoUaboration), PoS LAT2009, 079 (2009) farXiv: 0910. 3618 [hep-lat]]. 

[59] A. Bazavov et al. (MILC CoUaboration), PoS LAT2009, 077 (2009) |arXiv: 09 1 1.0472 [hep-lat]]. 

[60] K. Orginos and D. Toussaint (MILC Collaboration), Phys. Rev. D 59, 014501 l;i999) |arXiv:hep-lat/9805009| . 

[61] G. P. Lepage, Phys. Rev. D 59, 074502 (1999) arXiv:hep-lat/9809157 . 

[62] K. Orginos, D. Touss aint and R. L. Sugar XMILC Collaboration) , Phys. Rev. D 60, 054503 (1999) 

'arXiv : hep-lat/9903032' . 
[63] T. A. DeGrand, A. Hasenfratz and T. G. Kovacs, Phys. Rev. D 67, 054501 (2003) [arXiv : hep-lat/0211006 . 



